diff --git a/examples/t8code_2d_dgsem/elixir_three_equations_ellitptical_drop.jl b/examples/t8code_2d_dgsem/elixir_three_equations_ellitptical_drop.jl new file mode 100644 index 00000000000..e134c556d73 --- /dev/null +++ b/examples/t8code_2d_dgsem/elixir_three_equations_ellitptical_drop.jl @@ -0,0 +1,146 @@ +# The same setup as tree_2d_dgsem/elixir_advection_basic.jl +# to verify the StructuredMesh implementation against TreeMesh + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +gamma = 1.0 +k0 = 3.2e-5 +rho0 = 1000.0 +# G = 9.81 +G = 0.0 + +equations = ThreeEquations2D(gamma, k0, rho0, G) + +function smootherstep(left, right, x) + # Scale, and clamp x to 0..1 range. + x = clamp((x - left) / (right - left)) + + return x * x * x * (x * (6.0 * x - 15.0) + 10.0) +end + +@inline function clamp(x, lowerlimit = 0.0, upperlimit = 1.0) + if (x < lowerlimit) return lowerlimit end + if (x > upperlimit) return upperlimit end + return x +end + +function initial_condition(x, t, equations::ThreeEquations2D) + + r = Trixi.norm(x) + + width = 0.1 + radius = 1.0 + + s = smootherstep(radius + 0.5*width, radius - 0.5*width, r) + + rho = 1000.0 + v1 = -100.0 * x[1] * s + v2 = 100.0 * x[2] * s + alpha = clamp(s, 1e-3, 1.0) + + # # liquid domain + # if((x[1]^2 + x[2]^2) <= 1) + # rho = 1000.0 + # alpha = 1.0 - 10^(-3) + # v1 = -100.0 * x[1] + # v2 = 100.0 * x[2] + # else + # rho = 1000.0 + # v1 = 0.0 + # v2 = 0.0 + # alpha = 10^(-3) + # end + # phi = x[2] + phi = 0.0 + + # rho = 1000.0 + # v1 = 0.0 + # v2 = 0.0 + # alpha = 1.0 + + return prim2cons(SVector(rho, v1, v2, alpha, phi, 0.0), equations) +end + +# volume_flux = (flux_central, flux_nonconservative_ThreeEquations) +volume_flux = (Trixi.flux_entropy_cons_gamma_one_ThreeEquations, Trixi.flux_non_cons_entropy_cons_gamma_one_ThreeEquations) + +surface_flux = (Trixi.flux_entropy_cons_gamma_one_ThreeEquations, Trixi.flux_non_cons_entropy_cons_gamma_one_ThreeEquations) +# surface_flux = (flux_lax_friedrichs, Trixi.flux_non_cons_entropy_cons_gamma_one_ThreeEquations) + +polydeg = 3 + +basis = LobattoLegendreBasis(polydeg) + +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max=1.0, + alpha_min=0.001, + alpha_smooth=true, + variable=alpha_rho) + +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg=volume_flux, + volume_flux_fv=surface_flux) + +# volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux +# solver = DGSEM(basis, surface_flux = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-3.0, -3.0) # minimum coordinates (min(x), min(y)) +coordinates_max = ( 3.0, 3.0) # maximum coordinates (max(x), max(y)) + +trees_per_dimension = (1, 1) + +initial_refinement_level = 6 + +mesh = P4estMesh(trees_per_dimension, polydeg = polydeg , + coordinates_min = coordinates_min, coordinates_max = coordinates_max, + initial_refinement_level = initial_refinement_level) + +# A semidiscretization collects data structures and functions for the spatial discretization +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.0076) + +# Create ODE problem with time span from 0.0 to 1.0 +ode = semidiscretize(semi, tspan) + +# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup +# and resets the timers +summary_callback = SummaryCallback() + +# The SaveSolutionCallback allows to save the solution to a file in regular intervals +save_solution = SaveSolutionCallback(interval = 20, + solution_variables = cons2cons) + +# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results +analysis_callback = AnalysisCallback(semi, interval = 100) + +# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step +stepsize_callback = StepsizeCallback(cfl = 0.5) + +# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver +callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback, save_solution) + +############################################################################### +# run the simulation + +# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks +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, callback = callbacks); + +# Print the timer summary +summary_callback() + +# Finalize `T8codeMesh` to make sure MPI related objects in t8code are +# released before `MPI` finalizes. +!isinteractive() && finalize(mesh) diff --git a/src/equations/three_equations_2d.jl b/src/equations/three_equations_2d.jl index 7c47236f091..81861367915 100644 --- a/src/equations/three_equations_2d.jl +++ b/src/equations/three_equations_2d.jl @@ -260,6 +260,23 @@ return SVector(f1, f2, f3, f4, f5, f6) end + @inline function flux_nonconservative_ThreeEquations(u_ll, u_rr, orientation::Integer, + equations::ThreeEquations2D) + v1_ll = u_ll[2] / u_ll[1] + v2_ll = u_ll[3] / u_ll[1] + alpha_rr = u_rr[4] + + z = zero(eltype(u_ll)) + + if orientation == 1 + f = SVector(z, z, z, v1_ll * alpha_rr, z, z) + else + f = SVector(z, z, z, v2_ll * alpha_rr, z, z) + end + + return f + end + # Calculate 1D flux for a single point in the normal direction # Note, this directional vector is not normalized @inline function flux(u, normal_direction::AbstractVector, equations::ThreeEquations2D) @@ -280,40 +297,154 @@ return SVector(f1, f2, f3, f4, f5, f6) end - @inline function flux_nonconservative_ThreeEquations(u_ll, u_rr, orientation::Integer, + @inline function flux_nonconservative_ThreeEquations(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, equations::ThreeEquations2D) v1_ll = u_ll[2] / u_ll[1] v2_ll = u_ll[3] / u_ll[1] alpha_rr = u_rr[4] + v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] + z = zero(eltype(u_ll)) - if orientation == 1 - f = SVector(z, z, z, v1_ll * alpha_rr, z, z) - else - f = SVector(z, z, z, v2_ll * alpha_rr, z, z) - end + f = SVector(z, z, z, v_dot_n_ll * alpha_rr, z, z) return f end - @inline function flux_nonconservative_ThreeEquations(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, + @inline function flux_entropy_cons_ThreeEquations(u_ll, u_rr, + normal_direction::AbstractVector, equations::ThreeEquations2D) - v1_ll = u_ll[2] / u_ll[1] - v2_ll = u_ll[3] / u_ll[1] - alpha_rr = u_rr[4] + z = zero(eltype(u_ll)) - v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] + rho_ll, v1_ll, v2_ll, alpha_ll, phi_ll, dummy_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, alpha_rr, phi_rr, dummy_rr = cons2prim(u_rr, equations) + + rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma) + + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + alpha_avg = 0.5f0 * (alpha_ll + alpha_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] + + v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + + p_ll = pressure(u_ll, equations) + p_rr = pressure(u_rr, equations) + + p_avg = 0.5f0 * (p_ll + p_rr) + + f1 = alpha_avg * rho_mean * v_dot_n_avg + f2 = f1 * v1_avg + alpha_avg * p_avg * normal_direction[1] + f3 = f1 * v2_avg + alpha_avg * p_avg * normal_direction[2] + f4 = z + f5 = z + + return SVector(f1, f2, f3, f4, f5, z) + end + + @inline function flux_non_cons_entropy_cons_ThreeEquations(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ThreeEquations2D) z = zero(eltype(u_ll)) - f = SVector(z, z, z, v_dot_n_ll * alpha_rr, z, z) + rho_ll, v1_ll, v2_ll, alpha_ll, phi_ll, dummy_ll = cons2prim(u_ll, equations) + rho_rr, v1_rr, v2_rr, alpha_rr, phi_rr, dummy_rr = cons2prim(u_rr, equations) - return f + phi_jump = phi_rr - phi_ll + alpha_jump = alpha_rr - alpha_ll + + alpha_avg = 0.5f0 * (alpha_ll + alpha_rr) + + rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma) + # rho_mean = 0.5f0 * (rho_ll + rho_rr) + + v_dot_n_ll = v1_ll * normal_direction_average[1] + v2_ll * normal_direction_average[2] + v_dot_n_rr = v1_rr * normal_direction_average[1] + v2_rr * normal_direction_average[2] + + v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + + f1 = z + f2 = phi_jump * alpha_avg * rho_mean * normal_direction_average[1] + f3 = phi_jump * alpha_avg * rho_mean * normal_direction_average[2] + # f4 = alpha_jump * v_dot_n_avg + f4 = alpha_jump * v_dot_n_ll + f5 = z + + return SVector(f1, f2, f3, f4, f5, z) end + ## @inline function flux_entropy_cons_gamma_one_ThreeEquations(u_ll, u_rr, + ## normal_direction::AbstractVector, + ## equations::ThreeEquations2D) + ## z = zero(eltype(u_ll)) + + ## rho_ll, v1_ll, v2_ll, alpha_ll, phi_ll, dummy_ll = cons2prim(u_ll, equations) + ## rho_rr, v1_rr, v2_rr, alpha_rr, phi_rr, dummy_rr = cons2prim(u_rr, equations) + + ## rho_mean = ln_mean(rho_ll, rho_rr) + + ## v1_avg = 0.5f0 * (v1_ll + v1_rr) + ## v2_avg = 0.5f0 * (v2_ll + v2_rr) + + ## alpha_avg = 0.5f0 * (alpha_ll + alpha_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] + + ## v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + + ## p_ll = pressure(u_ll, equations) + ## p_rr = pressure(u_rr, equations) + + ## p_avg = 0.5f0 * (p_ll + p_rr) + + ## f1 = alpha_avg * rho_mean * v_dot_n_avg + ## f2 = f1 * v1_avg + alpha_avg * p_avg * normal_direction[1] + ## f3 = f1 * v2_avg + alpha_avg * p_avg * normal_direction[2] + ## f4 = z + ## f5 = z + + ## return SVector(f1, f2, f3, f4, f5, z) + ## end + + ## @inline function flux_non_cons_entropy_cons_gamma_one_ThreeEquations(u_ll, u_rr, + ## normal_direction_ll::AbstractVector, + ## normal_direction_average::AbstractVector, + ## equations::ThreeEquations2D) + ## z = zero(eltype(u_ll)) + + ## rho_ll, v1_ll, v2_ll, alpha_ll, phi_ll, dummy_ll = cons2prim(u_ll, equations) + ## rho_rr, v1_rr, v2_rr, alpha_rr, phi_rr, dummy_rr = cons2prim(u_rr, equations) + + ## phi_jump = phi_rr - phi_ll + ## alpha_jump = alpha_rr - alpha_ll + + ## alpha_avg = 0.5f0 * (alpha_ll + alpha_rr) + + ## rho_mean = ln_mean(rho_ll, rho_rr) + ## # rho_mean = 0.5f0 * (rho_ll + rho_rr) + + ## v_dot_n_ll = v1_ll * normal_direction_average[1] + v2_ll * normal_direction_average[2] + ## v_dot_n_rr = v1_rr * normal_direction_average[1] + v2_rr * normal_direction_average[2] + + ## v_dot_n_avg = 0.5f0 * (v_dot_n_ll + v_dot_n_rr) + + ## f1 = z + ## f2 = phi_jump * alpha_avg * rho_mean * normal_direction_average[1] + ## f3 = phi_jump * alpha_avg * rho_mean * normal_direction_average[2] + ## f4 = alpha_jump * v_dot_n_avg + ## f5 = z + + ## return SVector(f1, f2, f3, f4, f5, z) + ## end + @inline function flux_nonconservative_ThreeEquations_well(u_ll, u_rr, orientation::Integer, equations::ThreeEquations2D) v1_ll = u_ll[2] / u_ll[1] @@ -351,9 +482,9 @@ v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] # Non well-balanced force contribution. - force = -u_ll[1] * phi_rr + force = u_ll[1] * phi_rr - # Well-balanced force contribution. + # # Well-balanced force contribution. # force = -u_ll[1] * (equations.k0/equations.rho_0) * # exp(-equations.rho_0/equations.k0 * phi_ll) * # exp( equations.rho_0/equations.k0 * phi_rr) @@ -472,24 +603,17 @@ # Convert conservative variables to primitive @inline function cons2entropy(u, equations::ThreeEquations2D) - alpha_rho, alpha_rho_v1, alpha_rho_v2, alpha, phi, dummy = u - - rho = alpha_rho / alpha - v1 = alpha_rho_v1 / alpha_rho - v2 = alpha_rho_v2 / alpha_rho - - phi = max(1e-6, max(max_abs_speeds(u, equations)...) ) - - # vn = sqrt(v1*v1 + v2*v2) - # vn = log10(sqrt(alpha_rho_v1^2 + alpha_rho_v2^2)) - vn = log10(max(1e-6, sqrt(v1^2 + v2^2))) - - alpha_log = log10(max(1e-6, alpha)) + rho, v1, v2, alpha, phi, dummy = cons2prim(u, equations) p = pressure(u, equations) - # return SVector(alpha_rho, vn, log10(alpha), alpha, phi) - # return SVector(alpha_rho, vn, alpha_log, alpha, phi) - return SVector(alpha_rho, p, alpha_log, alpha, phi, dummy) + w1 = internal_energy_density(u, equations) - 0.5*(v1^2 + v2^2) + p/rho + phi + w2 = v1 + w3 = v2 + w4 = -p + w5 = alpha*rho + w6 = 0.0 + + return SVector(w1, w2, w3, w4, w5, w6) end # Convert primitive to conservative variables @@ -524,10 +648,37 @@ @inline function pressure(u, equations::ThreeEquations2D) alpha_rho, alpha_rho_v1, alpha_rho_v2, alpha, phi, dummy = u - # return alpha * equations.k0 * ((alpha_rho / alpha / equations.rho_0)^(equations.gamma) - 1) return equations.k0 * ((alpha_rho / alpha / equations.rho_0)^(equations.gamma) - 1) end + @inline function internal_energy_density(u, equations::ThreeEquations2D) + alpha_rho, alpha_rho_v1, alpha_rho_v2, alpha, phi, dummy = u + return alpha*equations.k0/alpha_rho * ((alpha_rho / alpha / equations.rho_0)^(equations.gamma)/(equations.gamma-1) + 1) + end + + @inline function internal_energy(u, equations::ThreeEquations2D) + alpha_rho, alpha_rho_v1, alpha_rho_v2, alpha, phi, dummy = u + return alpha_rho * internal_energy_density(u, equations) + end + + @inline function kinetic_energy(u, equations::ThreeEquations2D) + alpha_rho, alpha_rho_v1, alpha_rho_v2, alpha, phi, dummy = u + return 0.5*(alpha_rho_v1^2 + alpha_rho_v2^2)/alpha_rho + end + + @inline function potential_energy(u, equations::ThreeEquations2D) + alpha_rho, alpha_rho_v1, alpha_rho_v2, alpha, phi, dummy = u + return alpha_rho*phi + end + + @inline function total_energy(u, equations::ThreeEquations2D) + return kinetic_energy(u, equations) + internal_energy(u, equations) + potential_energy(u, equations) + end + + @inline function entropy(u, equations::ThreeEquations2D) + return total_energy(u, equations) + end + @inline function density_pressure(u, equations::ThreeEquations2D) alpha_rho, alpha_rho_v1, alpha_rho_v2, alpha, phi, dummy = u rho = alpha_rho / alpha