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.

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

Central SBP finite difference solver for UnstructuredMesh2D #1773

Merged
merged 7 commits into from
Dec 13, 2023
Merged
Changes from 1 commit
Commits
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
Prev Previous commit
Next Next commit
apply formatter to new and edited files
andrewwinters5000 committed Dec 13, 2023
commit 35ccce43a6276947125c6f450b64ec855976556e
43 changes: 22 additions & 21 deletions examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl
Original file line number Diff line number Diff line change
@@ -13,30 +13,31 @@ 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 )
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)
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())
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)
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)
@@ -45,7 +46,7 @@ mesh = UnstructuredMesh2D(mesh_file)
# create the semi discretization object

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

###############################################################################
# ODE solvers, callbacks etc.
@@ -56,13 +57,13 @@ ode = semidiscretize(semi, tspan)
summary_callback = SummaryCallback()

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

alive_callback = AliveCallback(analysis_interval=analysis_interval)
alive_callback = AliveCallback(analysis_interval = analysis_interval)

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

callbacks = CallbackSet(summary_callback, analysis_callback,
alive_callback, save_solution)
@@ -71,6 +72,6 @@ callbacks = CallbackSet(summary_callback, analysis_callback,
# 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)
sol = solve(ode, SSPRK43(), abstol = 1.0e-12, reltol = 1.0e-12,
save_everystep = false, callback = callbacks)
summary_callback() # print the timer summary
31 changes: 16 additions & 15 deletions examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl
Original file line number Diff line number Diff line change
@@ -14,27 +14,28 @@ 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)
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())
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)
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)
mesh = UnstructuredMesh2D(mesh_file, periodicity = true)

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

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

###############################################################################
# ODE solvers, callbacks etc.
@@ -45,20 +46,20 @@ ode = semidiscretize(semi, tspan)
summary_callback = SummaryCallback()

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

alive_callback = AliveCallback(analysis_interval=analysis_interval)
alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
save_initial_solution=true,
save_final_solution=true)
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)
sol = solve(ode, SSPRK43(), abstol = 1.0e-9, reltol = 1.0e-9,
save_everystep = false, callback = callbacks)
summary_callback() # print the timer summary
3 changes: 2 additions & 1 deletion src/solvers/dg.jl
Original file line number Diff line number Diff line change
@@ -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) && !(dg.volume_integral isa VolumeIntegralStrongForm)
if !(dg.volume_integral isa VolumeIntegralWeakForm) &&
!(dg.volume_integral isa VolumeIntegralStrongForm)
show(increment_indent(io), mime, dg.volume_integral)
end
summary_footer(io)
3 changes: 2 additions & 1 deletion src/solvers/dgsem_unstructured/containers_2d.jl
Original file line number Diff line number Diff line change
@@ -97,7 +97,8 @@ function init_elements!(elements::UnstructuredElementContainer2D, mesh, basis)
end

# initialize all the values in the container of a general element (either straight sided or curved)
function init_element!(elements, element, basis::LobattoLegendreBasis, corners_or_surface_curves)
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)

187 changes: 95 additions & 92 deletions src/solvers/fdsbp_unstructured/containers_2d.jl
Original file line number Diff line number Diff line change
@@ -6,116 +6,119 @@
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin
#! format: noindent

# initialize all the values in the container of a general FD block (either straight sided or curved)
# OBS! Requires the SBP derivative matrix in order to compute metric terms that are free-stream preserving
function init_element!(elements, element, basis::AbstractDerivativeOperator, corners_or_surface_curves)
function init_element!(elements, element, basis::AbstractDerivativeOperator,
corners_or_surface_curves)
calc_node_coordinates!(elements.node_coordinates, element, get_nodes(basis),
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, basis,
elements.node_coordinates)

calc_metric_terms!(elements.jacobian_matrix, element, basis, elements.node_coordinates)
calc_inverse_jacobian!(elements.inverse_jacobian, element, elements.jacobian_matrix)

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

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

calc_normal_directions!(elements.normal_directions, element, elements.jacobian_matrix)

return elements
return elements
end


# construct the metric terms for a FDSBP element "block". Directly use the derivative matrix
# applied to the node coordinates.
# TODO: FD; How to make this work for the upwind solver because basis has three available derivative matrices
function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeOperator, node_coordinates)

# storage format:
# jacobian_matrix[1,1,:,:,:] <- X_xi
# jacobian_matrix[1,2,:,:,:] <- X_eta
# jacobian_matrix[2,1,:,:,:] <- Y_xi
# jacobian_matrix[2,2,:,:,:] <- Y_eta

# Compute the xi derivatives by applying D on the left
# This is basically the same as
# jacobian_matrix[1, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[1, :, :, element]
# but uses only matrix-vector products instead of a matrix-matrix product.
for j in eachnode(D_SBP)
mul!(view(jacobian_matrix, 1, 1, :, j, element), D_SBP,
view(node_coordinates, 1, :, j, element))
end
# jacobian_matrix[2, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[2, :, :, element]
for j in eachnode(D_SBP)
mul!(view(jacobian_matrix, 2, 1, :, j, element), D_SBP,
view(node_coordinates, 2, :, j, element))
end

# Compute the eta derivatives by applying transpose of D on the right
# jacobian_matrix[1, 2, :, :, element] = node_coordinates[1, :, :, element] * Matrix(D_SBP)'
for i in eachnode(D_SBP)
mul!(view(jacobian_matrix, 1, 2, i, :, element), D_SBP,
view(node_coordinates, 1, i, :, element))
end
# jacobian_matrix[2, 2, :, :, element] = node_coordinates[2, :, :, element] * Matrix(D_SBP)'
for i in eachnode(D_SBP)
mul!(view(jacobian_matrix, 2, 2, i, :, element), D_SBP,
view(node_coordinates, 2, i, :, element))
end

return jacobian_matrix
function calc_metric_terms!(jacobian_matrix, element, D_SBP::AbstractDerivativeOperator,
node_coordinates)

# storage format:
# jacobian_matrix[1,1,:,:,:] <- X_xi
# jacobian_matrix[1,2,:,:,:] <- X_eta
# jacobian_matrix[2,1,:,:,:] <- Y_xi
# jacobian_matrix[2,2,:,:,:] <- Y_eta

# Compute the xi derivatives by applying D on the left
# This is basically the same as
# jacobian_matrix[1, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[1, :, :, element]
# but uses only matrix-vector products instead of a matrix-matrix product.
for j in eachnode(D_SBP)
mul!(view(jacobian_matrix, 1, 1, :, j, element), D_SBP,
view(node_coordinates, 1, :, j, element))
end
# jacobian_matrix[2, 1, :, :, element] = Matrix(D_SBP) * node_coordinates[2, :, :, element]
for j in eachnode(D_SBP)
mul!(view(jacobian_matrix, 2, 1, :, j, element), D_SBP,
view(node_coordinates, 2, :, j, element))
end

# Compute the eta derivatives by applying transpose of D on the right
# jacobian_matrix[1, 2, :, :, element] = node_coordinates[1, :, :, element] * Matrix(D_SBP)'
for i in eachnode(D_SBP)
mul!(view(jacobian_matrix, 1, 2, i, :, element), D_SBP,
view(node_coordinates, 1, i, :, element))
end
# jacobian_matrix[2, 2, :, :, element] = node_coordinates[2, :, :, element] * Matrix(D_SBP)'
for i in eachnode(D_SBP)
mul!(view(jacobian_matrix, 2, 2, i, :, element), D_SBP,
view(node_coordinates, 2, i, :, element))
end

return jacobian_matrix
end

# construct the normal direction vectors (but not actually normalized) for a curved sided FDSBP element "block"
# normalization occurs on the fly during the surface flux computation
# OBS! This assumes that the boundary points are included.
function calc_normal_directions!(normal_directions, element, jacobian_matrix)

# normal directions on the boundary for the left (local side 4) and right (local side 2)
N_ = size(jacobian_matrix, 3)
for j in 1:N_
# +x side or side 2 in the local indexing
X_xi = jacobian_matrix[1, 1, N_, j, element]
X_eta = jacobian_matrix[1, 2, N_, j, element]
Y_xi = jacobian_matrix[2, 1, N_, j, element]
Y_eta = jacobian_matrix[2, 2, N_, j, element]
Jtemp = X_xi * Y_eta - X_eta * Y_xi
normal_directions[1, j, 2, element] = sign(Jtemp) * ( Y_eta)
normal_directions[2, j, 2, element] = sign(Jtemp) * (-X_eta)

# -x side or side 4 in the local indexing
X_xi = jacobian_matrix[1, 1, 1, j, element]
X_eta = jacobian_matrix[1, 2, 1, j, element]
Y_xi = jacobian_matrix[2, 1, 1, j, element]
Y_eta = jacobian_matrix[2, 2, 1, j, element]
Jtemp = X_xi * Y_eta - X_eta * Y_xi
normal_directions[1, j, 4, element] = -sign(Jtemp) * ( Y_eta)
normal_directions[2, j, 4, element] = -sign(Jtemp) * (-X_eta)
end

# normal directions on the boundary for the top (local side 3) and bottom (local side 1)
N_ = size(jacobian_matrix, 4)
for i in 1:N_
# -y side or side 1 in the local indexing
X_xi = jacobian_matrix[1, 1, i, 1, element]
X_eta = jacobian_matrix[1, 2, i, 1, element]
Y_xi = jacobian_matrix[2, 1, i, 1, element]
Y_eta = jacobian_matrix[2, 2, i, 1, element]
Jtemp = X_xi * Y_eta - X_eta * Y_xi
normal_directions[1, i, 1, element] = -sign(Jtemp) * (-Y_xi)
normal_directions[2, i, 1, element] = -sign(Jtemp) * ( X_xi)

# +y side or side 3 in the local indexing
X_xi = jacobian_matrix[1, 1, i, N_, element]
X_eta = jacobian_matrix[1, 2, i, N_, element]
Y_xi = jacobian_matrix[2, 1, i, N_, element]
Y_eta = jacobian_matrix[2, 2, i, N_, element]
Jtemp = X_xi * Y_eta - X_eta * Y_xi
normal_directions[1, i, 3, element] = sign(Jtemp) * (-Y_xi)
normal_directions[2, i, 3, element] = sign(Jtemp) * ( X_xi)
end

return normal_directions
# normal directions on the boundary for the left (local side 4) and right (local side 2)
N_ = size(jacobian_matrix, 3)
andrewwinters5000 marked this conversation as resolved.
Show resolved Hide resolved
for j in 1:N_
# +x side or side 2 in the local indexing
X_xi = jacobian_matrix[1, 1, N_, j, element]
X_eta = jacobian_matrix[1, 2, N_, j, element]
Y_xi = jacobian_matrix[2, 1, N_, j, element]
Y_eta = jacobian_matrix[2, 2, N_, j, element]
Jtemp = X_xi * Y_eta - X_eta * Y_xi
normal_directions[1, j, 2, element] = sign(Jtemp) * (Y_eta)
normal_directions[2, j, 2, element] = sign(Jtemp) * (-X_eta)

# -x side or side 4 in the local indexing
X_xi = jacobian_matrix[1, 1, 1, j, element]
X_eta = jacobian_matrix[1, 2, 1, j, element]
Y_xi = jacobian_matrix[2, 1, 1, j, element]
Y_eta = jacobian_matrix[2, 2, 1, j, element]
Jtemp = X_xi * Y_eta - X_eta * Y_xi
normal_directions[1, j, 4, element] = -sign(Jtemp) * (Y_eta)
normal_directions[2, j, 4, element] = -sign(Jtemp) * (-X_eta)
end

# normal directions on the boundary for the top (local side 3) and bottom (local side 1)
N_ = size(jacobian_matrix, 4)
andrewwinters5000 marked this conversation as resolved.
Show resolved Hide resolved
for i in 1:N_
# -y side or side 1 in the local indexing
X_xi = jacobian_matrix[1, 1, i, 1, element]
X_eta = jacobian_matrix[1, 2, i, 1, element]
Y_xi = jacobian_matrix[2, 1, i, 1, element]
Y_eta = jacobian_matrix[2, 2, i, 1, element]
Jtemp = X_xi * Y_eta - X_eta * Y_xi
normal_directions[1, i, 1, element] = -sign(Jtemp) * (-Y_xi)
normal_directions[2, i, 1, element] = -sign(Jtemp) * (X_xi)

# +y side or side 3 in the local indexing
X_xi = jacobian_matrix[1, 1, i, N_, element]
X_eta = jacobian_matrix[1, 2, i, N_, element]
Y_xi = jacobian_matrix[2, 1, i, N_, element]
Y_eta = jacobian_matrix[2, 2, i, N_, element]
Jtemp = X_xi * Y_eta - X_eta * Y_xi
normal_directions[1, i, 3, element] = sign(Jtemp) * (-Y_xi)
normal_directions[2, i, 3, element] = sign(Jtemp) * (X_xi)
end

return normal_directions
end


end # @muladd
end # @muladd
Loading