Skip to content

Commit

Permalink
apply formatter to new and edited files
Browse files Browse the repository at this point in the history
  • Loading branch information
andrewwinters5000 committed Dec 13, 2023
1 parent 5d2c04a commit 35ccce4
Show file tree
Hide file tree
Showing 8 changed files with 289 additions and 282 deletions.
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
Expand Up @@ -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)
Expand All @@ -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.
Expand All @@ -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)
Expand All @@ -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
Expand Up @@ -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.
Expand All @@ -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
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) && !(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)
Expand Down
3 changes: 2 additions & 1 deletion src/solvers/dgsem_unstructured/containers_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
187 changes: 95 additions & 92 deletions src/solvers/fdsbp_unstructured/containers_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
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
end


end # @muladd
end # @muladd
Loading

0 comments on commit 35ccce4

Please sign in to comment.