Skip to content

Commit

Permalink
add option to subtract a reference state in SaveSolutionCallback
Browse files Browse the repository at this point in the history
  • Loading branch information
@tristanmontoya
tristanmontoya committed Oct 24, 2024
1 parent 4a81e00 commit f2fbe1c
Show file tree
Hide file tree
Showing 3 changed files with 102 additions and 37 deletions.
58 changes: 53 additions & 5 deletions examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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))

Expand All @@ -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)
output_directory = "out_bubble",
solution_variables = cons2prim,
reference_solution = background_state)

stepsize_callback = StepsizeCallback(cfl = 1.0)

Expand All @@ -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
Expand Down
44 changes: 41 additions & 3 deletions src/callbacks_step/save_solution.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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})
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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

Expand Down
37 changes: 8 additions & 29 deletions src/callbacks_step/save_solution_dg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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,
Expand Down

0 comments on commit f2fbe1c

Please sign in to comment.