diff --git a/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl b/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl new file mode 100644 index 0000000000..c181203e7a --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl @@ -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 diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl new file mode 100644 index 0000000000..7ada50c0c6 --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl @@ -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 diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl new file mode 100644 index 0000000000..edcd221bf5 --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl @@ -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 diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 9e5ebc7f9b..9b61df62cc 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -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. """ @@ -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) @@ -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). diff --git a/src/solvers/dgsem_unstructured/containers_2d.jl b/src/solvers/dgsem_unstructured/containers_2d.jl index 13eeaeabff..f51dd09801 100644 --- a/src/solvers/dgsem_unstructured/containers_2d.jl +++ b/src/solvers/dgsem_unstructured/containers_2d.jl @@ -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) @@ -84,24 +84,25 @@ 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) @@ -109,7 +110,7 @@ function init_element!(elements, element, nodes, corners_or_surface_curves) 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 diff --git a/src/solvers/fdsbp_tree/fdsbp_2d.jl b/src/solvers/fdsbp_tree/fdsbp_2d.jl index beff605629..09d18cecd7 100644 --- a/src/solvers/fdsbp_tree/fdsbp_2d.jl +++ b/src/solvers/fdsbp_tree/fdsbp_2d.jl @@ -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), diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl new file mode 100644 index 0000000000..3857c2d8a2 --- /dev/null +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -0,0 +1,124 @@ +# !!! warning "Experimental implementation (curvilinear FDSBP)" +# This is an experimental feature and may change in future releases. + +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# 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) + 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_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, + elements.jacobian_matrix) + + 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 +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, 4) + 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, 3) + 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 diff --git a/src/solvers/fdsbp_unstructured/fdsbp.jl b/src/solvers/fdsbp_unstructured/fdsbp.jl new file mode 100644 index 0000000000..dee9776abb --- /dev/null +++ b/src/solvers/fdsbp_unstructured/fdsbp.jl @@ -0,0 +1,14 @@ +# !!! warning "Experimental implementation (curvilinear FDSBP)" +# This is an experimental feature and may change in future releases. + +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +# dimension specific curvilinear implementations and data structures +include("containers_2d.jl") +include("fdsbp_2d.jl") +end # @muladd diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl new file mode 100644 index 0000000000..b459f4c42c --- /dev/null +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -0,0 +1,219 @@ +# !!! warning "Experimental implementation (curvilinear FDSBP)" +# This is an experimental feature and may change in future releases. + +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +# 2D unstructured cache +function create_cache(mesh::UnstructuredMesh2D, equations, dg::FDSBP, RealT, uEltype) + elements = init_elements(mesh, equations, dg.basis, RealT, uEltype) + + interfaces = init_interfaces(mesh, elements) + + boundaries = init_boundaries(mesh, elements) + + cache = (; elements, interfaces, boundaries) + + # Add specialized parts of the cache required to for efficient flux computations + cache = (; cache..., + create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...) + + return cache +end + +# TODO: FD; Upwind versions of surface / volume integral + +# 2D volume integral contributions for `VolumeIntegralStrongForm` +# OBS! This is the standard (not de-aliased) form of the volume integral. +# So it is not provably stable for variable coefficients due to the the metric terms. +@inline function calc_volume_integral!(du, u, + mesh::UnstructuredMesh2D, + nonconservative_terms::False, equations, + volume_integral::VolumeIntegralStrongForm, + dg::FDSBP, cache) + D = dg.basis # SBP derivative operator + @unpack f_threaded = cache + @unpack contravariant_vectors = cache.elements + + # SBP operators from SummationByPartsOperators.jl implement the basic interface + # of matrix-vector multiplication. Thus, we pass an "array of structures", + # packing all variables per node in an `SVector`. + if nvariables(equations) == 1 + # `reinterpret(reshape, ...)` removes the leading dimension only if more + # than one variable is used. + u_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(u)}, u), + nnodes(dg), nnodes(dg), nelements(dg, cache)) + du_vectors = reshape(reinterpret(SVector{nvariables(equations), eltype(du)}, + du), + nnodes(dg), nnodes(dg), nelements(dg, cache)) + else + u_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(u)}, u) + du_vectors = reinterpret(reshape, SVector{nvariables(equations), eltype(du)}, + du) + end + + # Use the tensor product structure to compute the discrete derivatives of + # the contravariant fluxes line-by-line and add them to `du` for each element. + @threaded for element in eachelement(dg, cache) + f_element = f_threaded[Threads.threadid()] + u_element = view(u_vectors, :, :, element) + + # x direction + for j in eachnode(dg) + for i in eachnode(dg) + Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element) + f_element[i, j] = flux(u_element[i, j], Ja1, equations) + end + mul!(view(du_vectors, :, j, element), D, view(f_element, :, j), + one(eltype(du)), one(eltype(du))) + end + + # y direction + for i in eachnode(dg) + for j in eachnode(dg) + Ja2 = get_contravariant_vector(2, contravariant_vectors, i, j, element) + f_element[i, j] = flux(u_element[i, j], Ja2, equations) + end + mul!(view(du_vectors, i, :, element), D, view(f_element, i, :), + one(eltype(du)), one(eltype(du))) + end + end + + return nothing +end + +# Note! The local side numbering for the unstructured quadrilateral element implementation differs +# from the structured TreeMesh or StructuredMesh local side numbering: +# +# TreeMesh/StructuredMesh sides versus UnstructuredMesh sides +# 4 3 +# ----------------- ----------------- +# | | | | +# | ^ eta | | ^ eta | +# 1 | | | 2 4 | | | 2 +# | | | | | | +# | ---> xi | | ---> xi | +# ----------------- ----------------- +# 3 1 +# Therefore, we require a different surface integral routine here despite their similar structure. +# Also, the normal directions are already outward pointing for `UnstructuredMesh2D` so all the +# surface contributions are added. +function calc_surface_integral!(du, u, mesh::UnstructuredMesh2D, + equations, surface_integral::SurfaceIntegralStrongForm, + dg::DG, cache) + inv_weight_left = inv(left_boundary_weight(dg.basis)) + inv_weight_right = inv(right_boundary_weight(dg.basis)) + @unpack normal_directions, surface_flux_values = cache.elements + + @threaded for element in eachelement(dg, cache) + for l in eachnode(dg) + # surface at -x + u_node = get_node_vars(u, equations, dg, 1, l, element) + # compute internal flux in normal direction on side 4 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 4, + element) + f_node = flux(u_node, outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 4, element) + multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), + equations, dg, 1, l, element) + + # surface at +x + u_node = get_node_vars(u, equations, dg, nnodes(dg), l, element) + # compute internal flux in normal direction on side 2 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 2, + element) + f_node = flux(u_node, outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 2, element) + multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), + equations, dg, nnodes(dg), l, element) + + # surface at -y + u_node = get_node_vars(u, equations, dg, l, 1, element) + # compute internal flux in normal direction on side 1 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 1, + element) + f_node = flux(u_node, outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 1, element) + multiply_add_to_node_vars!(du, inv_weight_left, (f_num - f_node), + equations, dg, l, 1, element) + + # surface at +y + u_node = get_node_vars(u, equations, dg, l, nnodes(dg), element) + # compute internal flux in normal direction on side 3 + outward_direction = get_node_coords(normal_directions, equations, dg, l, 3, + element) + f_node = flux(u_node, outward_direction, equations) + f_num = get_node_vars(surface_flux_values, equations, dg, l, 3, element) + multiply_add_to_node_vars!(du, inv_weight_right, (f_num - f_node), + equations, dg, l, nnodes(dg), element) + end + end + + return nothing +end + +# AnalysisCallback +function integrate_via_indices(func::Func, u, + mesh::UnstructuredMesh2D, equations, + dg::FDSBP, cache, args...; normalize = true) where {Func} + # TODO: FD. This is rather inefficient right now and allocates... + weights = diag(SummationByPartsOperators.mass_matrix(dg.basis)) + + # Initialize integral with zeros of the right shape + integral = zero(func(u, 1, 1, 1, equations, dg, args...)) + total_volume = zero(real(mesh)) + + # Use quadrature to numerically integrate over entire domain + for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, element])) + integral += volume_jacobian * weights[i] * weights[j] * + func(u, i, j, element, equations, dg, args...) + total_volume += volume_jacobian * weights[i] * weights[j] + end + end + + # Normalize with total volume + if normalize + integral = integral / total_volume + end + + return integral +end + +function calc_error_norms(func, u, t, analyzer, + mesh::UnstructuredMesh2D, equations, initial_condition, + dg::FDSBP, cache, cache_analysis) + # TODO: FD. This is rather inefficient right now and allocates... + weights = diag(SummationByPartsOperators.mass_matrix(dg.basis)) + @unpack node_coordinates, inverse_jacobian = cache.elements + + # Set up data structures + l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1, 1), equations)) + linf_error = copy(l2_error) + total_volume = zero(real(mesh)) + + # Iterate over all elements for error calculations + for element in eachelement(dg, cache) + for j in eachnode(analyzer), i in eachnode(analyzer) + volume_jacobian = abs(inv(cache.elements.inverse_jacobian[i, j, element])) + u_exact = initial_condition(get_node_coords(node_coordinates, equations, dg, + i, j, element), t, equations) + diff = func(u_exact, equations) - + func(get_node_vars(u, equations, dg, i, j, element), equations) + l2_error += diff .^ 2 * (weights[i] * weights[j] * volume_jacobian) + linf_error = @. max(linf_error, abs(diff)) + total_volume += weights[i] * weights[j] * volume_jacobian + end + end + + # For L2 error, divide by total volume + l2_error = @. sqrt(l2_error / total_volume) + + return l2_error, linf_error +end +end # @muladd diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 5341d86a7d..139b423ead 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -664,6 +664,67 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +# TODO: FD; for now put the unstructured tests for the 2D FDSBP here. +@trixi_testset "FDSBP (central): elixir_advection_basic.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"), + "elixir_advection_basic.jl"), + l2=[0.0001105211407319266], + linf=[0.0004199363734466166]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "FDSBP (central): elixir_euler_source_terms.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"), + "elixir_euler_source_terms.jl"), + l2=[8.155544666380138e-5, + 0.0001477863788446318, + 0.00014778637884460072, + 0.00045584189984542687], + linf=[0.0002670775876922882, + 0.0005683064706873964, + 0.0005683064706762941, + 0.0017770812025146299], + tspan=(0.0, 0.05)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "FDSBP (central): elixir_euler_free_stream.jl" begin + @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"), + "elixir_euler_free_stream.jl"), + l2=[5.4329175009362306e-14, + 1.0066867437607972e-13, + 6.889210012578449e-14, + 1.568290814572709e-13], + linf=[5.963762816918461e-10, + 5.08869890669672e-11, + 1.1581377523661729e-10, + 4.61017890529547e-11], + tspan=(0.0, 0.1), + atol=1.0e-11) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end end # Clean up afterwards: delete Trixi.jl output directory