diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl index 38b9386e94e..cd1986dcced 100644 --- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl +++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl @@ -26,6 +26,22 @@ struct WarmBubbleSetup end end +@inline function cons2potentialtemperature(u, equations::CompressibleEulerEquations2D) + p_0 = 1.0f5 # Pascals + c_p = 1004.0f0 + c_v = 717.0f0 + R = c_p - c_v + rho, v1, v2, p = cons2prim(u, equations) + exner_pressure = (p / p_0)^(R / c_p) + T = p / (rho * R) + potential_temperature = T / exner_pressure + return SVector(rho, v1, v2, potential_temperature) +end + +function Trixi.varnames(::typeof(cons2potentialtemperature), ::CompressibleEulerEquations2D) + return ("rho", "v1", "v2", "Potential temperature") +end + # Initial condition function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) @unpack g, c_p, c_v = setup @@ -73,8 +89,36 @@ end return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) end +# Background reference state +function hydrostatic_background(x, t, equations::CompressibleEulerEquations2D) + g = 9.81 + c_p = 1004.0f0 + c_v = 717.0f0 + potential_temperature = 300.0 + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] + + # pressure + p_0 = 100_000.0 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # temperature + T = potential_temperature * exner + + # density + rho = p / (R * T) + + v1 = 20.0 + v2 = 0.0 + E = c_v * T + 0.5 * (v1^2 + v2^2) + return SVector(rho, rho * v1, rho * v2, rho * E) +end + ############################################################################### # semidiscretization of the compressible Euler equations + warm_bubble_setup = WarmBubbleSetup() equations = CompressibleEulerEquations2D(warm_bubble_setup.gamma) @@ -92,14 +136,14 @@ basis = LobattoLegendreBasis(polydeg) surface_flux = FluxLMARS(340.0) volume_flux = flux_kennedy_gruber -volume_integral = VolumeIntegralFluxDifferencing(volume_flux) - +#volume_integral = VolumeIntegralFluxDifferencing(volume_flux) +volume_integral = VolumeIntegralWeakForm() solver = DGSEM(basis, surface_flux, volume_integral) coordinates_min = (0.0, 0.0) coordinates_max = (20_000.0, 10_000.0) -cells_per_dimension = (64, 32) +cells_per_dimension = (200, 100) mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, periodicity = (true, false)) @@ -123,11 +167,14 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, alive_callback = AliveCallback(analysis_interval = analysis_interval) +background_state = compute_reference_state(hydrostatic_background, semi, cons2prim) + save_solution = SaveSolutionCallback(interval = analysis_interval, save_initial_solution = true, save_final_solution = true, output_directory = "out", - solution_variables = cons2prim) + solution_variables = cons2prim, + reference_solution = background_state) stepsize_callback = StepsizeCallback(cfl = 1.0) @@ -139,6 +186,7 @@ callbacks = CallbackSet(summary_callback, ############################################################################### # run the simulation + sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), maxiters = 1.0e7, dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback diff --git a/src/Trixi.jl b/src/Trixi.jl index 0cedec782df..5d5d223d5b5 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -268,6 +268,8 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback, AnalysisSurfaceIntegral, DragCoefficientPressure, LiftCoefficientPressure, DragCoefficientShearStress, LiftCoefficientShearStress +export compute_reference_state + export load_mesh, load_time, load_timestep, load_timestep!, load_dt, load_adaptive_time_integrator! diff --git a/src/callbacks_step/save_solution.jl b/src/callbacks_step/save_solution.jl index 3f9abef7f33..1853cf2ce69 100644 --- a/src/callbacks_step/save_solution.jl +++ b/src/callbacks_step/save_solution.jl @@ -21,12 +21,14 @@ at a single point to a set of solution variables. The first parameter passed to `solution_variables` will be the set of conservative variables and the second parameter is the equation struct. """ -mutable struct SaveSolutionCallback{IntervalType, SolutionVariablesType} +mutable struct SaveSolutionCallback{IntervalType, SolutionVariablesType, + ReferenceSolutionType} interval_or_dt::IntervalType save_initial_solution::Bool save_final_solution::Bool output_directory::String solution_variables::SolutionVariablesType + reference_solution::ReferenceSolutionType end function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SaveSolutionCallback}) @@ -96,7 +98,8 @@ function SaveSolutionCallback(; interval::Integer = 0, save_initial_solution = true, save_final_solution = true, output_directory = "out", - solution_variables = cons2prim) + solution_variables = cons2prim, + reference_solution = nothing) if !isnothing(dt) && interval > 0 throw(ArgumentError("You can either set the number of steps between output (using `interval`) or the time between outputs (using `dt`) but not both simultaneously")) end @@ -110,7 +113,8 @@ function SaveSolutionCallback(; interval::Integer = 0, solution_callback = SaveSolutionCallback(interval_or_dt, save_initial_solution, save_final_solution, - output_directory, solution_variables) + output_directory, solution_variables, + reference_solution) # Expected most frequent behavior comes first if isnothing(dt) @@ -243,6 +247,40 @@ end node_variables; system = system) end +# Convenience function to convert variables +function convert_variables(u, equations, solution_variables) + if solution_variables === cons2cons + data = u + n_vars = nvariables(equations) + else + # Reinterpret the solution array as an array of conservative variables, + # compute the solution variables via broadcasting, and reinterpret the + # result as a plain array of floating point numbers + data = Array(reinterpret(eltype(u), + solution_variables.(reinterpret(SVector{nvariables(equations), + eltype(u)}, u), + Ref(equations)))) + + # Find out variable count by looking at output from `solution_variables` function + n_vars = size(data, 1) + end + + return data, n_vars +end + +# Compute a reference state to be subtracted from the solution at each write to file +function compute_reference_state(func, semi, solution_variables) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + + reference_state = compute_coefficients(func, 0.0, semi) + reference_state_wrapped = copy(Trixi.wrap_array_native(reference_state, + mesh, equations, solver, + cache)) + + data, _ = convert_variables(reference_state_wrapped, equations, solution_variables) + return data +end + # TODO: Taal refactor, move save_mesh_file? # function save_mesh_file(mesh::TreeMesh, output_directory, timestep=-1) in io/io.jl diff --git a/src/callbacks_step/save_solution_dg.jl b/src/callbacks_step/save_solution_dg.jl index 57cdb3ef8a2..9c556487a9a 100644 --- a/src/callbacks_step/save_solution_dg.jl +++ b/src/callbacks_step/save_solution_dg.jl @@ -15,7 +15,7 @@ function save_solution_file(u, time, dt, timestep, element_variables = Dict{Symbol, Any}(), node_variables = Dict{Symbol, Any}(); system = "") - @unpack output_directory, solution_variables = solution_callback + @unpack output_directory, solution_variables, reference_solution = solution_callback # Filename based on current time step if isempty(system) @@ -26,20 +26,13 @@ function save_solution_file(u, time, dt, timestep, end # Convert to different set of variables if requested - if solution_variables === cons2cons - data = u - n_vars = nvariables(equations) + converted_variables, n_vars = convert_variables(u, equations, solution_variables) + + # Subtract reference solution + if !isnothing(reference_solution) + data = converted_variables - reference_solution else - # Reinterpret the solution array as an array of conservative variables, - # compute the solution variables via broadcasting, and reinterpret the - # result as a plain array of floating point numbers - data = Array(reinterpret(eltype(u), - solution_variables.(reinterpret(SVector{nvariables(equations), - eltype(u)}, u), - Ref(equations)))) - - # Find out variable count by looking at output from `solution_variables` function - n_vars = size(data, 1) + data = converted_variables end # Open file (clobber existing content) @@ -108,21 +101,7 @@ function save_solution_file(u, time, dt, timestep, end # Convert to different set of variables if requested - if solution_variables === cons2cons - data = u - n_vars = nvariables(equations) - else - # Reinterpret the solution array as an array of conservative variables, - # compute the solution variables via broadcasting, and reinterpret the - # result as a plain array of floating point numbers - data = Array(reinterpret(eltype(u), - solution_variables.(reinterpret(SVector{nvariables(equations), - eltype(u)}, u), - Ref(equations)))) - - # Find out variable count by looking at output from `solution_variables` function - n_vars = size(data, 1) - end + data, n_vars = convert_variables(u, equations, solution_variables) if HDF5.has_parallel() save_solution_file_parallel(data, time, dt, timestep, n_vars, mesh, equations,