From 40f6547a2fc227ce4c10d15da9a03790b430f0bb Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 5 Mar 2024 21:40:42 +0100 Subject: [PATCH] Apply suggestions from code review Co-authored-by: Hendrik Ranocha --- .../elixir_euler_free_stream_upwind.jl | 8 ++++++-- .../elixir_euler_source_terms_upwind.jl | 10 +++++++--- src/solvers/fdsbp_unstructured/fdsbp_2d.jl | 11 +++++++---- 3 files changed, 20 insertions(+), 9 deletions(-) diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl index 4f391f81f89..2a1956f9d10 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_free_stream_upwind.jl @@ -1,3 +1,5 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# This is an experimental feature and may change in future releases. using OrdinaryDiffEq using Trixi @@ -20,6 +22,7 @@ boundary_conditions = Dict(:outerCircle => boundary_condition_free_stream, ############################################################################### # Get the Upwind FDSBP approximation space +# TODO: FDSBP # Note, one must set `xmin=-1` and `xmax=1` due to the reuse # of interpolation routines from `calc_node_coordinates!` to create # the physical coordinates in the mappings. @@ -37,8 +40,9 @@ solver = FDSBP(D_upw, ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) -# Mesh with second order boundary polynomials requires an upwind SBP operator -# with (at least) 4th order boundary closure to guarantee the approximation is free-stream +# Mesh with second-order boundary polynomials requires an upwind SBP operator +# with (at least) 4th order boundary closure to guarantee the approximation is +# free-stream preserving mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/ec9a345f09199ebe471d35d5c1e4e08f/raw/15975943d8642e42f8292235314b6f1b30aa860d/mesh_inner_outer_boundaries.mesh", joinpath(@__DIR__, "mesh_inner_outer_boundaries.mesh")) diff --git a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl index 8a4cc70e764..3bb2c4f1cc9 100644 --- a/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl +++ b/examples/unstructured_2d_fdsbp/elixir_euler_source_terms_upwind.jl @@ -1,3 +1,5 @@ +# !!! warning "Experimental implementation (upwind SBP)" +# This is an experimental feature and may change in future releases. using OrdinaryDiffEq using Trixi @@ -21,6 +23,7 @@ boundary_conditions = Dict(:Top => boundary_condition_eoc, ############################################################################### # Get the Upwind FDSBP approximation space +# TODO: FDSBP # Note, one must set `xmin=-1` and `xmax=1` due to the reuse # of interpolation routines from `calc_node_coordinates!` to create # the physical coordinates in the mappings. @@ -38,15 +41,16 @@ solver = FDSBP(D_upw, ############################################################################### # Get the curved quad mesh from a file (downloads the file if not available locally) -# Mesh with first order boundary polynomials requires an upwind SBP operator -# with (at least) 2nd order boundary closure to guarantee the approximation is free-stream +# Mesh with first-order boundary polynomials requires an upwind SBP operator +# with (at least) 2nd order boundary closure to guarantee the approximation is +# free-stream preserving mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/a4f4743008bf3233957a9ea6ac7a62e0/raw/8b36cc6649153fe0a5723b200368a210a1d74eaf/mesh_refined_box.mesh", joinpath(@__DIR__, "mesh_refined_box.mesh")) mesh = UnstructuredMesh2D(mesh_file) ############################################################################### -# create the semi discretization object +# create the semidiscretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_term, diff --git a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl index 948c020a1fb..f0aca6ec690 100644 --- a/src/solvers/fdsbp_unstructured/fdsbp_2d.jl +++ b/src/solvers/fdsbp_unstructured/fdsbp_2d.jl @@ -125,16 +125,19 @@ function calc_volume_integral!(du, u, # the fluxes line-by-line and add them to `du` for each element. @threaded for element in eachelement(dg, cache) # f_minus_plus_element wraps the storage provided by f_minus_element and - # f_plus_element such that we can use a single plain broadcasting below. - # f_minus_element and f_plus_element are updated in broadcasting calls - # of the form `@. f_minus_plus_element = ...`. + # f_plus_element such that we can use a single assignment below. + # f_minus_element and f_plus_element are updated whenever we update + # `f_minus_plus_element[i, j] = ...` below. f_minus_plus_element = f_minus_plus_threaded[Threads.threadid()] f_minus_element = f_minus_threaded[Threads.threadid()] f_plus_element = f_plus_threaded[Threads.threadid()] u_element = view(u_vectors, :, :, element) # x direction - # @. f_minus_plus_element = splitting(u_element, 1, equations) + # We use flux vector splittings in the directions of the contravariant + # basis vectors. Thus, we do not use a broadcasting operation like + # @. f_minus_plus_element = splitting(u_element, 1, equations) + # in the Cartesian case but loop over all nodes. for j in eachnode(dg), i in eachnode(dg) # contravariant vectors computed with central D matrix Ja1 = get_contravariant_vector(1, contravariant_vectors, i, j, element)