Skip to content

Commit

Permalink
Merge branch 'main' into simplessprk_dt
Browse files Browse the repository at this point in the history
  • Loading branch information
amrueda authored Dec 13, 2023
2 parents b8d1a2a + 14796ea commit 4844491
Show file tree
Hide file tree
Showing 11 changed files with 644 additions and 12 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Trixi"
uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
version = "0.6.4-pre"
version = "0.6.5-pre"

[deps]
CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
Expand Down
69 changes: 69 additions & 0 deletions examples/unstructured_2d_fdsbp/elixir_advection_basic.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@

using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advection_velocity = (0.2, -0.7)
equations = LinearScalarAdvectionEquation2D(advection_velocity)

###############################################################################
# Get the FDSBP approximation operator

D_SBP = derivative_operator(SummationByPartsOperators.MattssonAlmquistVanDerWeide2018Accurate(),
derivative_order = 1, accuracy_order = 4,
xmin = -1.0, xmax = 1.0, N = 15)
solver = FDSBP(D_SBP,
surface_integral = SurfaceIntegralStrongForm(flux_lax_friedrichs),
volume_integral = VolumeIntegralStrongForm())

###############################################################################
# Get the curved quad mesh from a file (downloads the file if not available locally)

default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")
isfile(default_mesh_file) ||
download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh",
default_mesh_file)
mesh_file = default_mesh_file

mesh = UnstructuredMesh2D(mesh_file, periodicity = true)

###############################################################################
# create the semidiscretization object

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
solver)

###############################################################################
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.0
ode = semidiscretize(semi, (0.0, 1.0))

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval = 100)

# The SaveSolutionCallback allows to save the solution to a file in regular intervals
save_solution = SaveSolutionCallback(interval = 100,
solution_variables = cons2prim)

# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl = 1.6)

# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_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, callback = callbacks);
summary_callback() # print the timer summary
77 changes: 77 additions & 0 deletions examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,77 @@

using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

# Free-stream initial condition
initial_condition = initial_condition_constant

# Boundary conditions for free-stream testing
boundary_condition_free_stream = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = Dict(:Body => boundary_condition_free_stream,
:Button1 => boundary_condition_free_stream,
:Button2 => boundary_condition_free_stream,
:Eye1 => boundary_condition_free_stream,
:Eye2 => boundary_condition_free_stream,
:Smile => boundary_condition_free_stream,
:Bowtie => boundary_condition_free_stream)

###############################################################################
# Get the FDSBP approximation space

D_SBP = derivative_operator(SummationByPartsOperators.MattssonAlmquistVanDerWeide2018Accurate(),
derivative_order = 1, accuracy_order = 4,
xmin = -1.0, xmax = 1.0, N = 12)
solver = FDSBP(D_SBP,
surface_integral = SurfaceIntegralStrongForm(flux_hll),
volume_integral = VolumeIntegralStrongForm())

###############################################################################
# Get the curved quad mesh from a file (downloads the file if not available locally)

default_mesh_file = joinpath(@__DIR__, "mesh_gingerbread_man.mesh")
isfile(default_mesh_file) ||
download("https://gist.githubusercontent.com/andrewwinters5000/2c6440b5f8a57db131061ad7aa78ee2b/raw/1f89fdf2c874ff678c78afb6fe8dc784bdfd421f/mesh_gingerbread_man.mesh",
default_mesh_file)
mesh_file = default_mesh_file

mesh = UnstructuredMesh2D(mesh_file)

###############################################################################
# create the semi discretization object

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 5.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 100,
save_initial_solution = true,
save_final_solution = true)

callbacks = CallbackSet(summary_callback, analysis_callback,
alive_callback, save_solution)

###############################################################################
# run the simulation

# set small tolerances for the free-stream preservation test
sol = solve(ode, SSPRK43(), abstol = 1.0e-12, reltol = 1.0e-12,
save_everystep = false, callback = callbacks)
summary_callback() # print the timer summary
65 changes: 65 additions & 0 deletions examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,65 @@

using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

initial_condition = initial_condition_convergence_test

###############################################################################
# Get the FDSBP approximation operator

D_SBP = derivative_operator(SummationByPartsOperators.MattssonNordström2004(),
derivative_order = 1, accuracy_order = 4,
xmin = -1.0, xmax = 1.0, N = 10)
solver = FDSBP(D_SBP,
surface_integral = SurfaceIntegralStrongForm(flux_lax_friedrichs),
volume_integral = VolumeIntegralStrongForm())

###############################################################################
# Get the curved quad mesh from a file (downloads the file if not available locally)

default_mesh_file = joinpath(@__DIR__, "mesh_periodic_square_with_twist.mesh")
isfile(default_mesh_file) ||
download("https://gist.githubusercontent.com/andrewwinters5000/12ce661d7c354c3d94c74b964b0f1c96/raw/8275b9a60c6e7ebbdea5fc4b4f091c47af3d5273/mesh_periodic_square_with_twist.mesh",
default_mesh_file)
mesh_file = default_mesh_file

mesh = UnstructuredMesh2D(mesh_file, periodicity = true)

###############################################################################
# create the semi discretization object

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
source_terms = source_terms_convergence_test)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 100,
save_initial_solution = true,
save_final_solution = true)

callbacks = CallbackSet(summary_callback, analysis_callback,
alive_callback, save_solution)

###############################################################################
# run the simulation

sol = solve(ode, SSPRK43(), abstol = 1.0e-9, reltol = 1.0e-9,
save_everystep = false, callback = callbacks)
summary_callback() # print the timer summary
8 changes: 5 additions & 3 deletions src/solvers/dg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -41,8 +41,8 @@ standard textbooks.
Applications
[doi: 10.1007/978-0-387-72067-8](https://doi.org/10.1007/978-0-387-72067-8)
`VolumeIntegralWeakForm()` is only implemented for conserved terms as
non-conservative terms should always be discretized in conjunction with a flux-splitting scheme,
`VolumeIntegralWeakForm()` is only implemented for conserved terms as
non-conservative terms should always be discretized in conjunction with a flux-splitting scheme,
see [`VolumeIntegralFluxDifferencing`](@ref).
This treatment is required to achieve, e.g., entropy-stability or well-balancedness.
"""
Expand Down Expand Up @@ -415,7 +415,8 @@ function Base.show(io::IO, mime::MIME"text/plain", dg::DG)
summary_line(io, "surface integral", dg.surface_integral |> typeof |> nameof)
show(increment_indent(io), mime, dg.surface_integral)
summary_line(io, "volume integral", dg.volume_integral |> typeof |> nameof)
if !(dg.volume_integral isa VolumeIntegralWeakForm)
if !(dg.volume_integral isa VolumeIntegralWeakForm) &&
!(dg.volume_integral isa VolumeIntegralStrongForm)
show(increment_indent(io), mime, dg.volume_integral)
end
summary_footer(io)
Expand Down Expand Up @@ -598,6 +599,7 @@ include("dgsem/dgsem.jl")
# and boundary conditions weakly. Thus, these methods can re-use a lot of
# functionality implemented for DGSEM.
include("fdsbp_tree/fdsbp.jl")
include("fdsbp_unstructured/fdsbp.jl")

function allocate_coefficients(mesh::AbstractMesh, equations, dg::DG, cache)
# We must allocate a `Vector` in order to be able to `resize!` it (AMR).
Expand Down
15 changes: 8 additions & 7 deletions src/solvers/dgsem_unstructured/containers_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ end
eachelement(elements::UnstructuredElementContainer2D)
Return an iterator over the indices that specify the location in relevant data structures
for the elements in `elements`.
for the elements in `elements`.
In particular, not the elements themselves are returned.
"""
@inline function eachelement(elements::UnstructuredElementContainer2D)
Expand Down Expand Up @@ -84,32 +84,33 @@ function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis)
# loop through elements and call the correct constructor based on whether the element is curved
for element in eachelement(elements)
if mesh.element_is_curved[element]
init_element!(elements, element, basis.nodes,
init_element!(elements, element, basis,
view(mesh.surface_curves, :, element))
else # straight sided element
for i in 1:4, j in 1:2
# pull the (x,y) values of these corners out of the global corners array
four_corners[i, j] = mesh.corners[j, mesh.element_node_ids[i, element]]
end
init_element!(elements, element, basis.nodes, four_corners)
init_element!(elements, element, basis, four_corners)
end
end
end

# initialize all the values in the container of a general element (either straight sided or curved)
function init_element!(elements, element, nodes, corners_or_surface_curves)
calc_node_coordinates!(elements.node_coordinates, element, nodes,
function init_element!(elements, element, basis::LobattoLegendreBasis,
corners_or_surface_curves)
calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis),
corners_or_surface_curves)

calc_metric_terms!(elements.jacobian_matrix, element, nodes,
calc_metric_terms!(elements.jacobian_matrix, element, get_nodes(basis),
corners_or_surface_curves)

calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix)

calc_contravariant_vectors!(elements.contravariant_vectors, element,
elements.jacobian_matrix)

calc_normal_directions!(elements.normal_directions, element, nodes,
calc_normal_directions!(elements.normal_directions, element, get_nodes(basis),
corners_or_surface_curves)

return elements
Expand Down
2 changes: 1 addition & 1 deletion src/solvers/fdsbp_tree/fdsbp_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
#! format: noindent

# 2D caches
function create_cache(mesh::TreeMesh{2}, equations,
function create_cache(mesh::Union{TreeMesh{2}, UnstructuredMesh2D}, equations,
volume_integral::VolumeIntegralStrongForm, dg, uEltype)
prototype = Array{SVector{nvariables(equations), uEltype}, ndims(mesh)}(undef,
ntuple(_ -> nnodes(dg),
Expand Down
Loading

0 comments on commit 4844491

Please sign in to comment.