Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: Add option to subtract a reference state in SaveSolutionCallback #2131

Draft
wants to merge 2 commits into
base: main
Choose a base branch
from
Draft

WIP: Add option to subtract a reference state in SaveSolutionCallback #2131

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f2fbe1c
add option to subtract a reference state in SaveSolutionCallback
tristanmontoya Oct 24, 2024
3b86205
export compute_reference_state
tristanmontoya Oct 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 52 additions & 4 deletions examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
View file
Open in desktop
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)
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
2 changes: 2 additions & 0 deletions src/Trixi.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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!

Expand Down
44 changes: 41 additions & 3 deletions src/callbacks_step/save_solution.jl
View file
Open in desktop
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
Open in desktop
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
Loading