From 806374b0c98c71c4d89d00186e8d5fb411a054f9 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Dec 2023 11:55:52 +0100 Subject: [PATCH 1/6] containers and kernels for FDSBP solver on UnstructuredMesh2D --- src/solvers/dg.jl | 7 +- .../dgsem_unstructured/containers_2d.jl | 14 +- src/solvers/fdsbp_tree/fdsbp_2d.jl | 2 +- .../fdsbp_unstructured/containers_2d.jl | 121 ++++++++++ src/solvers/fdsbp_unstructured/fdsbp.jl | 16 ++ src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 217 ++++++++++++++++++ 6 files changed, 366 insertions(+), 11 deletions(-) create mode 100644 src/solvers/fdsbp_unstructured/containers_2d.jl create mode 100644 src/solvers/fdsbp_unstructured/fdsbp.jl create mode 100644 src/solvers/fdsbp_unstructured/fdsbp_2d.jl diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 9e5ebc7f9b5..58d15102fa2 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,7 @@ 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 +598,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 13eeaeabffb..7995547cf18 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,24 @@ 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 +109,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 beff605629a..09d18cecd75 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 00000000000..7296f2ea73f --- /dev/null +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -0,0 +1,121 @@ +# !!! 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 + +# 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, 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 \ No newline at end of file diff --git a/src/solvers/fdsbp_unstructured/fdsbp.jl b/src/solvers/fdsbp_unstructured/fdsbp.jl new file mode 100644 index 00000000000..8645d08de82 --- /dev/null +++ b/src/solvers/fdsbp_unstructured/fdsbp.jl @@ -0,0 +1,16 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# 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 + + +# dimension specific curvilinear implementations and data structures +include("containers_2d.jl") +include("fdsbp_2d.jl") + + +end # @muladd \ No newline at end of file diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl new file mode 100644 index 00000000000..6384e7dbf23 --- /dev/null +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -0,0 +1,217 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# 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 + + +# 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 \ No newline at end of file From 5d2c04a6b916dc2f07825315af67d3461eb0bda9 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Dec 2023 11:56:27 +0100 Subject: [PATCH 2/6] add elixirs and corresponding tests --- .../elixir_euler_free_stream.jl | 76 +++++++++++++++++++ .../elixir_euler_source_terms.jl | 64 ++++++++++++++++ test/test_unstructured_2d.jl | 46 +++++++++++ 3 files changed, 186 insertions(+) create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl create mode 100644 examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl 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 00000000000..7b17353c73a --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl @@ -0,0 +1,76 @@ + +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 00000000000..023e436d573 --- /dev/null +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl @@ -0,0 +1,64 @@ + +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/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 5341d86a7d1..bc6be6b98e9 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -664,6 +664,52 @@ 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_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 From 35ccce43a6276947125c6f450b64ec855976556e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Dec 2023 12:12:51 +0100 Subject: [PATCH 3/6] apply formatter to new and edited files --- .../elixir_euler_free_stream.jl | 43 +-- .../elixir_euler_source_terms.jl | 31 +- src/solvers/dg.jl | 3 +- .../dgsem_unstructured/containers_2d.jl | 3 +- .../fdsbp_unstructured/containers_2d.jl | 187 +++++------ src/solvers/fdsbp_unstructured/fdsbp.jl | 8 +- src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 294 +++++++++--------- test/test_unstructured_2d.jl | 2 +- 8 files changed, 289 insertions(+), 282 deletions(-) diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl index 7b17353c73a..7ada50c0c65 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream.jl @@ -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 diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl index 023e436d573..edcd221bf59 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms.jl @@ -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,13 +46,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) @@ -59,6 +60,6 @@ callbacks = CallbackSet(summary_callback, analysis_callback, ############################################################################### # 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 diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index 58d15102fa2..9b61df62cc3 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -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) diff --git a/src/solvers/dgsem_unstructured/containers_2d.jl b/src/solvers/dgsem_unstructured/containers_2d.jl index 7995547cf18..f51dd09801b 100644 --- a/src/solvers/dgsem_unstructured/containers_2d.jl +++ b/src/solvers/dgsem_unstructured/containers_2d.jl @@ -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) diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl index 7296f2ea73f..7196747384d 100644 --- a/src/solvers/fdsbp_unstructured/containers_2d.jl +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -6,63 +6,68 @@ # 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" @@ -70,52 +75,50 @@ end # 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 \ No newline at end of file +end # @muladd diff --git a/src/solvers/fdsbp_unstructured/fdsbp.jl b/src/solvers/fdsbp_unstructured/fdsbp.jl index 8645d08de82..dee9776abb7 100644 --- a/src/solvers/fdsbp_unstructured/fdsbp.jl +++ b/src/solvers/fdsbp_unstructured/fdsbp.jl @@ -1,4 +1,4 @@ -# !!! warning "Experimental implementation (upwind SBP)" +# !!! 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). @@ -6,11 +6,9 @@ # 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 \ No newline at end of file +end # @muladd diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl index 6384e7dbf23..b459f4c42cc 100644 --- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -1,4 +1,4 @@ -# !!! warning "Experimental implementation (upwind SBP)" +# !!! 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). @@ -6,26 +6,25 @@ # 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) - elements = init_elements(mesh, equations, dg.basis, RealT, uEltype) - - interfaces = init_interfaces(mesh, elements) + interfaces = init_interfaces(mesh, elements) - boundaries = init_boundaries(mesh, elements) + boundaries = init_boundaries(mesh, elements) - cache = (; elements, interfaces, boundaries) + 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)...) + # 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 + return cache end - # TODO: FD; Upwind versions of surface / volume integral # 2D volume integral contributions for `VolumeIntegralStrongForm` @@ -36,56 +35,57 @@ end 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))) + 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 - # 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))) + # 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 - end - return nothing + return nothing end - # Note! The local side numbering for the unstructured quadrilateral element implementation differs # from the structured TreeMesh or StructuredMesh local side numbering: # @@ -105,113 +105,115 @@ end 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) + 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 - end - return nothing + 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] + 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 - end - # Normalize with total volume - if normalize - integral = integral / total_volume - end + # Normalize with total volume + if normalize + integral = integral / total_volume + end - return integral + 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 + # 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 - end - # For L2 error, divide by total volume - l2_error = @. sqrt(l2_error / total_volume) + # For L2 error, divide by total volume + l2_error = @. sqrt(l2_error / total_volume) - return l2_error, linf_error + return l2_error, linf_error end - -end # @muladd \ No newline at end of file +end # @muladd diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index bc6be6b98e9..9c22c05a332 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -700,7 +700,7 @@ end 1.1581377523661729e-10, 4.61017890529547e-11], tspan=(0.0, 0.1), - atol = 1.0e-11) + atol=1.0e-11) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 37b4bd8ea6d0511d288b1e878a69c2f3b1f000c5 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Dec 2023 13:59:59 +0100 Subject: [PATCH 4/6] add advection equation test to up coverage --- .../elixir_advection_basic.jl | 69 +++++++++++++++++++ test/test_unstructured_2d.jl | 15 ++++ 2 files changed, 84 insertions(+) create mode 100644 examples/unstructured_2d_fdsbp/elixir_advection_basic.jl 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 00000000000..d1944a6320d --- /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 semi discretization 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/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 9c22c05a332..139b423ead1 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -666,6 +666,21 @@ 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"), From 316555db60fbdbf9d1debb445e69f14dce3bddbd Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Dec 2023 15:37:07 +0100 Subject: [PATCH 5/6] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- examples/unstructured_2d_fdsbp/elixir_advection_basic.jl | 4 ++-- src/solvers/fdsbp_unstructured/containers_2d.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl b/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl index d1944a6320d..c181203e7a4 100644 --- a/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl +++ b/examples/unstructured_2d_fdsbp/elixir_advection_basic.jl @@ -31,7 +31,7 @@ mesh_file = default_mesh_file mesh = UnstructuredMesh2D(mesh_file, periodicity = true) ############################################################################### -# create the semi discretization object +# create the semidiscretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver) @@ -40,7 +40,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergen # ODE solvers, callbacks etc. # Create ODE problem with time span from 0.0 to 1.0 -ode = semidiscretize(semi, (0.0, 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 diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl index 7196747384d..b02507841db 100644 --- a/src/solvers/fdsbp_unstructured/containers_2d.jl +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -76,7 +76,7 @@ end 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) + 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] @@ -98,7 +98,7 @@ function calc_normal_directions!(normal_directions, element, jacobian_matrix) end # normal directions on the boundary for the top (local side 3) and bottom (local side 1) - N_ = size(jacobian_matrix, 4) + 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] From 4cae6f39e6db7674ecd92bfd47eeced421c2f35e Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Wed, 13 Dec 2023 15:38:36 +0100 Subject: [PATCH 6/6] update variable name to N --- .../fdsbp_unstructured/containers_2d.jl | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/solvers/fdsbp_unstructured/containers_2d.jl b/src/solvers/fdsbp_unstructured/containers_2d.jl index b02507841db..3857c2d8a20 100644 --- a/src/solvers/fdsbp_unstructured/containers_2d.jl +++ b/src/solvers/fdsbp_unstructured/containers_2d.jl @@ -77,12 +77,12 @@ 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_ + 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] + 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) @@ -99,7 +99,7 @@ function calc_normal_directions!(normal_directions, element, jacobian_matrix) # 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_ + 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] @@ -110,10 +110,10 @@ function calc_normal_directions!(normal_directions, element, jacobian_matrix) 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] + 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)