From ba70d67c3089a6a21786dcc22a35a4cdb428bf5d Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 08:54:43 -0500 Subject: [PATCH 01/35] minor change for consistency --- src/callbacks_step/glm_speed_dg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_step/glm_speed_dg.jl b/src/callbacks_step/glm_speed_dg.jl index dce64d5d04..eef01ed047 100644 --- a/src/callbacks_step/glm_speed_dg.jl +++ b/src/callbacks_step/glm_speed_dg.jl @@ -23,7 +23,7 @@ function calc_dt_for_cleaning_speed(cfl::Real, mesh, # Compute time step for GLM linear advection equation with c_h=1 for a DGMulti discretization. # Copies implementation behavior of `calc_dt_for_cleaning_speed` for DGSEM discretizations. - max_scaled_speed_for_c_h = (1 / minimum(md.J)) * ndims(equations) + max_scaled_speed_for_c_h = inv(minimum(md.J)) * ndims(equations) # This mimics `max_dt` for `TreeMesh`, except that `nnodes(dg)` is replaced by # `polydeg+1`. This is because `nnodes(dg)` returns the total number of From 51c2a5e3d61293a5af34e744a9a463d1c1207285 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 08:54:56 -0500 Subject: [PATCH 02/35] formatting --- src/solvers/dgsem_tree/dg_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index d7f1463fde..445a8082ce 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -702,7 +702,7 @@ end function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:Any,4}, t, boundary_condition, nonconservative_terms::True, equations, - surface_integral ,dg::DG, cache, + surface_integral, dg::DG, cache, direction, first_boundary, last_boundary) surface_flux, nonconservative_flux = surface_integral.surface_flux @unpack u, neighbor_ids, neighbor_sides, node_coordinates, orientations = cache.boundaries From 059c0e86725c471d95a742d9ab1fa4b7e0c15ee6 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 09:21:29 -0500 Subject: [PATCH 03/35] add nonconservative terms to DGMulti `calc_boundary_flux!` --- src/solvers/dgmulti/dg.jl | 15 ++++++++++----- src/solvers/dgmulti/flux_differencing.jl | 4 +++- .../dgmulti/flux_differencing_gauss_sbp.jl | 1 + 3 files changed, 14 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index f9e30f8f87..e1363219e0 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -428,23 +428,27 @@ end # do nothing for periodic (default) boundary conditions calc_boundary_flux!(cache, t, boundary_conditions::BoundaryConditionPeriodic, + have_nonconservative_terms, mesh, equations, dg::DGMulti) = nothing # "lispy tuple programming" instead of for loop for type stability -function calc_boundary_flux!(cache, t, boundary_conditions, mesh, equations, dg::DGMulti) +function calc_boundary_flux!(cache, t, boundary_conditions, have_nonconservative_terms, + mesh, equations, dg::DGMulti) # peel off first boundary condition calc_single_boundary_flux!(cache, t, first(boundary_conditions), first(keys(boundary_conditions)), - mesh, equations, dg) + have_nonconservative_terms, mesh, equations, dg) # recurse on the remainder of the boundary conditions - calc_boundary_flux!(cache, t, Base.tail(boundary_conditions), mesh, equations, dg) + calc_boundary_flux!(cache, t, Base.tail(boundary_conditions), + have_nonconservative_terms, mesh, equations, dg) end # terminate recursion calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple{(),Tuple{}}, - mesh, equations, dg::DGMulti) = nothing + have_nonconservative_terms, mesh, equations, dg::DGMulti) = nothing function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, + have_nonconservative_terms::False, mesh, equations, dg::DGMulti{NDIMS}) where {NDIMS} rd = dg.basis @@ -568,7 +572,8 @@ function rhs!(du, u, t, mesh, equations, have_nonconservative_terms(equations), equations, dg) @trixi_timeit timer() "boundary flux" calc_boundary_flux!( - cache, t, boundary_conditions, mesh, equations, dg) + cache, t, boundary_conditions, + have_nonconservative_terms(equations), mesh, equations, dg) @trixi_timeit timer() "surface integral" calc_surface_integral!( du, u, dg.surface_integral, mesh, equations, dg, cache) diff --git a/src/solvers/dgmulti/flux_differencing.jl b/src/solvers/dgmulti/flux_differencing.jl index 97905f1d0b..b541941388 100644 --- a/src/solvers/dgmulti/flux_differencing.jl +++ b/src/solvers/dgmulti/flux_differencing.jl @@ -595,6 +595,7 @@ function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions: equations, dg) @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, + have_nonconservative_terms(equations), mesh, equations, dg) @trixi_timeit timer() "surface integral" calc_surface_integral!(du, u, dg.surface_integral, @@ -630,7 +631,8 @@ function rhs!(du, u, t, mesh, equations, have_nonconservative_terms(equations), equations, dg) @trixi_timeit timer() "boundary flux" calc_boundary_flux!( - cache, t, boundary_conditions, mesh, equations, dg) + cache, t, boundary_conditions, have_nonconservative_terms(equations), + mesh, equations, dg) @trixi_timeit timer() "surface integral" calc_surface_integral!( du, u, dg.surface_integral, mesh, equations, dg, cache) diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index dfafe4ff98..b4542c527f 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl @@ -508,6 +508,7 @@ function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions: equations, dg) @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, + have_nonconservative_terms(equations), mesh, equations, dg) # `du` is stored at Gauss nodes here From 7d820b9e8a10494df483a7a54196e42f8b8a41f3 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 09:39:23 -0500 Subject: [PATCH 04/35] add noncon boundary flux --- src/solvers/dgmulti/dg.jl | 48 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index e1363219e0..1e0fd7ce9f 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -489,6 +489,54 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, # However, we don't have to re-reshape, since cache.flux_face_values still retains its original shape. end +function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, + have_nonconservative_terms::True, + mesh, equations, dg::DGMulti{NDIMS}) where {NDIMS} + + rd = dg.basis + md = mesh.md + (; u_face_values, flux_face_values) = cache + (; xyzf, nxyzJ, Jf) = md + surface_flux, nonconservative_flux = surface_integral.surface_flux + + # reshape face/normal arrays to have size = (num_points_on_face, num_faces_total). + # mesh.boundary_faces indexes into the columns of these face-reshaped arrays. + num_pts_per_face = rd.Nfq ÷ rd.Nfaces + num_faces_total = rd.Nfaces * md.num_elements + + # This function was originally defined as + # `reshape_by_face(u) = reshape(view(u, :), num_pts_per_face, num_faces_total)`. + # This results in allocations due to https://github.com/JuliaLang/julia/issues/36313. + # To avoid allocations, we use Tim Holy's suggestion: + # https://github.com/JuliaLang/julia/issues/36313#issuecomment-782336300. + reshape_by_face(u) = Base.ReshapedArray(u, (num_pts_per_face, num_faces_total), ()) + + u_face_values = reshape_by_face(u_face_values) + flux_face_values = reshape_by_face(flux_face_values) + Jf = reshape_by_face(Jf) + nxyzJ, xyzf = reshape_by_face.(nxyzJ), reshape_by_face.(xyzf) # broadcast over nxyzJ::NTuple{NDIMS,Matrix} + + # loop through boundary faces, which correspond to columns of reshaped u_face_values, ... + for f in mesh.boundary_faces[boundary_key] + for i in Base.OneTo(num_pts_per_face) + face_normal = SVector{NDIMS}(getindex.(nxyzJ, i, f)) / Jf[i,f] + face_coordinates = SVector{NDIMS}(getindex.(xyzf, i, f)) + + # compute conservative and non-conservative parts separately + flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, + surface_flux, equations) * Jf[i,f] + + flux_face_values[i,f] = flux_at_face_node + + 0.5 * boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, + nonconservative_flux, equations) * Jf[i,f] + + end + end + + # Note: modifying the values of the reshaped array modifies the values of cache.flux_face_values. + # However, we don't have to re-reshape, since cache.flux_face_values still retains its original shape. +end + # inverts Jacobian and scales by -1.0 function invert_jacobian!(du, mesh::DGMultiMesh, equations, dg::DGMulti, cache; scaling=-1) From e2d14329c0abae5d5d003d4b912bba76ae468962 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 09:42:11 -0500 Subject: [PATCH 05/35] fix dropped dg.surface_flux --- src/solvers/dgmulti/dg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 1e0fd7ce9f..083bd8ead1 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -497,7 +497,7 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, md = mesh.md (; u_face_values, flux_face_values) = cache (; xyzf, nxyzJ, Jf) = md - surface_flux, nonconservative_flux = surface_integral.surface_flux + surface_flux, nonconservative_flux = dg.surface_integral.surface_flux # reshape face/normal arrays to have size = (num_points_on_face, num_faces_total). # mesh.boundary_faces indexes into the columns of these face-reshaped arrays. From cafb3f174401b13158dd89806b780cec41a16afd Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 11:17:36 -0500 Subject: [PATCH 06/35] formatting --- src/solvers/dgmulti/dg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 083bd8ead1..8ea49b9b4a 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -524,7 +524,7 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, # compute conservative and non-conservative parts separately flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, - surface_flux, equations) * Jf[i,f] + surface_flux, equations) * Jf[i,f] flux_face_values[i,f] = flux_at_face_node + 0.5 * boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, From 8987e9a8ac5c1fb1c3391e2b5491f2aa343ea63f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 12:45:08 -0500 Subject: [PATCH 07/35] clean up noncons BCs --- src/solvers/dgmulti/dg.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 8ea49b9b4a..6c8dff4368 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -523,12 +523,13 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, face_coordinates = SVector{NDIMS}(getindex.(xyzf, i, f)) # compute conservative and non-conservative parts separately - flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, - surface_flux, equations) * Jf[i,f] + cons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, + surface_flux, equations) - flux_face_values[i,f] = flux_at_face_node + - 0.5 * boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, - nonconservative_flux, equations) * Jf[i,f] + noncons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, + nonconservative_flux, equations) + + flux_face_values[i,f] = (cons_flux_at_face_node + 0.5 * noncons_flux_at_face_node) * Jf[i,f] end end From 92ac9be2e7ee2c57b1127a98efb127be8e2dea5b Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 12:57:25 -0500 Subject: [PATCH 08/35] adding specialization of nonconservative Powell flux for BCs --- src/equations/ideal_glm_mhd_2d.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index c19273737e..f51e8f652d 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -238,6 +238,11 @@ terms. return f end +# TODO: is there a way to do this without manually dispatching for every non-conservative flux? We use this for the evaluation of boundary conditions. +# if only one normal is specified, use that for `normal_direction_average` +@inline flux_nonconservative_powell(u_ll, u_rr, normal_direction::AbstractVector, equations::IdealGlmMhdEquations2D) = + flux_nonconservative_powell(u_ll, u_rr, normal_direction, normal_direction, equations) + @inline function flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll::AbstractVector, normal_direction_average::AbstractVector, From ba9801dd74ff3a952d0c826be4ee0797fb61ca74 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 12:58:45 -0500 Subject: [PATCH 09/35] fix BoundaryConditionDoNothing for nonconservative terms --- src/basic_types.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/basic_types.jl b/src/basic_types.jl index ea7d542c5d..b7bcede275 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -77,19 +77,22 @@ struct BoundaryConditionDoNothing end @inline function (::BoundaryConditionDoNothing)( u_inner, orientation_or_normal_direction, direction::Integer, x, t, surface_flux, equations) return flux(u_inner, orientation_or_normal_direction, equations) + # TODO: should this be switched to `surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)` for consistency? end # This version can be called by hyperbolic solvers on unstructured, curved meshes @inline function (::BoundaryConditionDoNothing)( u_inner, outward_direction::AbstractVector, x, t, surface_flux, equations) - return flux(u_inner, outward_direction, equations) + + # this should reduce to `flux(u_inner, outward_direction, equations)` for a consistent symmetric flux. + return surface_flux(u_inner, u_inner, outward_direction, equations) end # This version can be called by parabolic solvers @inline function (::BoundaryConditionDoNothing)(inner_flux_or_state, other_args...) return inner_flux_or_state end - + """ boundary_condition_do_nothing = Trixi.BoundaryConditionDoNothing() From 22899b40ec6cdaf417d2b207398acd61bd22c180 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 12:58:50 -0500 Subject: [PATCH 10/35] add elixir --- .../dgmulti_2d/elixir_mhd_reflective_BCs.jl | 99 +++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl new file mode 100644 index 0000000000..5cfb3c6d01 --- /dev/null +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -0,0 +1,99 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the compressible ideal GLM-MHD equations +equations = IdealGlmMhdEquations2D(1.4) + +function initial_condition_perturbation(x, t, equations::IdealGlmMhdEquations2D) + # pressure perturbation in a vertically magnetized field on the domain [-1, 1]^2 + + r2 = (x[1] + 0.25)^2 + (x[2] + 0.25)^2 + + rho = 1.0 + v1 = 0.0 + v2 = 0.0 + v3 = 0.0 + p = 1 + 0.5 * exp(-100 * r2) + + # the pressure and magnetic field are chosen to be strongly + # magnetized, such that p / ||B||^2 ≈ 0.01. + B1 = 0.0 + B2 = 40.0 / sqrt(4.0 * pi) + B3 = 0.0 + + psi = 0.0 + return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations) +end +initial_condition = initial_condition_perturbation + +surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell) +volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) + +solver = DGMulti(polydeg=3, element_type = Quad(), approximation_type = SBP(), + surface_integral = SurfaceIntegralWeakForm(surface_flux), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +x_neg(x, tol=50*eps()) = abs(x[1] + 1) < tol +x_pos(x, tol=50*eps()) = abs(x[1] - 1) < tol +y_neg(x, tol=50*eps()) = abs(x[2] + 1) < tol +y_pos(x, tol=50*eps()) = abs(x[2] - 1) < tol +is_on_boundary = Dict(:x_neg => x_neg, :x_pos => x_pos, :y_neg => y_neg, :y_pos => y_pos) + +cells_per_dimension = (16, 16) +mesh = DGMultiMesh(solver, cells_per_dimension; periodicity=(false, false), is_on_boundary) + +function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector, x, t, + surface_flux_function, equations::IdealGlmMhdEquations2D) + + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + norm_ = norm(normal_direction) + normal = normal_direction / norm_ + + # compute the primitive variables + rho, v1, v2, v3, p, B1, B2, B3, psi = cons2prim(u_inner, equations) + + v_normal = dot(normal, SVector(v1, v2)) + u_mirror = prim2cons(SVector(rho, v1 - 2 * v_normal * normal[1], + v2 - 2 * v_normal * normal[2], + v3, p, B1, B2, B3, psi), equations) + + return surface_flux_function(u_inner, u_mirror, normal, equations) * norm_ +end + +boundary_conditions = (; x_neg=boundary_condition_velocity_slip_wall, + x_pos=boundary_condition_velocity_slip_wall, + y_neg=boundary_condition_do_nothing, + y_pos=BoundaryConditionDirichlet(initial_condition)) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; + boundary_conditions=boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, .075) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +alive_callback = AliveCallback(alive_interval=10) + +cfl = 0.5 +stepsize_callback = StepsizeCallback(cfl=cfl) +glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) + +callbacks = CallbackSet(summary_callback, + alive_callback, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), + dt=1e-5, # 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 From 24a4f1e0534ea9fbee07076fa87a377c2ce7c641 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 13:04:20 -0500 Subject: [PATCH 11/35] add test --- test/test_dgmulti_2d.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index d2eaeea57d..9730486af4 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -292,6 +292,14 @@ isdir(outdir) && rm(outdir, recursive=true) ) end + @trixi_testset "elixir_mhd_reflective_BCs.jl (Quad)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_reflective_BCs.jl"), + cells_per_dimension = 4, element_type = Quad(), approximation_type = GaussSBP(), + l2 = [0.0036019562526885497, 0.0017340971255526598, 0.008375221676923516, 0.0, 0.028596802653952254, 0.0018573697892179616, 0.002080779894047767, 0.0, 5.3012597624420295e-5], + linf = [0.01692598382370969, 0.00936965952976002, 0.04145170727841342, 0.0, 0.11569901084179435, 0.009849648257862952, 0.011417088537134745, 0.0, 0.0002992621756978016] + ) + end + @trixi_testset "elixir_shallowwater_source_terms.jl (Quad, SBP)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), cells_per_dimension = 8, element_type = Quad(), approximation_type = SBP(), From 0596c329007227b439e6d7816f64d3aea3e66b37 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 13:07:28 -0500 Subject: [PATCH 12/35] comment --- examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index 5cfb3c6d01..4c31aec0d4 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -44,6 +44,8 @@ is_on_boundary = Dict(:x_neg => x_neg, :x_pos => x_pos, :y_neg => y_neg, :y_pos cells_per_dimension = (16, 16) mesh = DGMultiMesh(solver, cells_per_dimension; periodicity=(false, false), is_on_boundary) +# Create a "reflective-like" boundary condition by mirroring the velocity but leaving the magnetic field alone. +# Note that this boundary condition is probably not entropy stable. function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector, x, t, surface_flux_function, equations::IdealGlmMhdEquations2D) From e05b7311728ec7de0099a41614b204b4ab0af209 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 14:32:52 -0500 Subject: [PATCH 13/35] importing norm --- examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index 4c31aec0d4..dfbf4c7929 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -46,6 +46,7 @@ mesh = DGMultiMesh(solver, cells_per_dimension; periodicity=(false, false), is_o # Create a "reflective-like" boundary condition by mirroring the velocity but leaving the magnetic field alone. # Note that this boundary condition is probably not entropy stable. +using LinearAlgebra: norm function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector, x, t, surface_flux_function, equations::IdealGlmMhdEquations2D) From f7131e47979ef8779f7cf94e3d69d5ce7ad81365 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 15:08:06 -0500 Subject: [PATCH 14/35] import dot as well --- examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index dfbf4c7929..205c33eec6 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -46,7 +46,7 @@ mesh = DGMultiMesh(solver, cells_per_dimension; periodicity=(false, false), is_o # Create a "reflective-like" boundary condition by mirroring the velocity but leaving the magnetic field alone. # Note that this boundary condition is probably not entropy stable. -using LinearAlgebra: norm +using LinearAlgebra: norm, dot function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector, x, t, surface_flux_function, equations::IdealGlmMhdEquations2D) From aefa0cd93bfc4bc103043d7653d29a637a4f24cf Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 16:15:35 -0500 Subject: [PATCH 15/35] adding forgotten analysis callback --- examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index 205c33eec6..cd623109ed 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -81,6 +81,8 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) alive_callback = AliveCallback(alive_interval=10) cfl = 0.5 @@ -88,6 +90,7 @@ stepsize_callback = StepsizeCallback(cfl=cfl) glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, stepsize_callback, glm_speed_callback) From 271399a732cad89b973ec59a729d8703eb1f1dc8 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Mon, 8 May 2023 16:16:07 -0500 Subject: [PATCH 16/35] Update src/solvers/dgmulti/dg.jl Co-authored-by: Michael Schlottke-Lakemper --- src/solvers/dgmulti/dg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 6c8dff4368..85f31cb275 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -501,8 +501,8 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, # reshape face/normal arrays to have size = (num_points_on_face, num_faces_total). # mesh.boundary_faces indexes into the columns of these face-reshaped arrays. - num_pts_per_face = rd.Nfq ÷ rd.Nfaces - num_faces_total = rd.Nfaces * md.num_elements + num_pts_per_face = rd.Nfq ÷ StartUpDG.num_faces(rd.element_type) + num_faces_total = StartUpDG.num_faces(rd.element_type) * md.num_elements # This function was originally defined as # `reshape_by_face(u) = reshape(view(u, :), num_pts_per_face, num_faces_total)`. From 383eb242bf1ea596042cef133bc70e6f70da0edb Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 16:17:29 -0500 Subject: [PATCH 17/35] remove some name-based type instabilities --- src/solvers/dgmulti/dg.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 85f31cb275..f6c387d3db 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -495,8 +495,6 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, rd = dg.basis md = mesh.md - (; u_face_values, flux_face_values) = cache - (; xyzf, nxyzJ, Jf) = md surface_flux, nonconservative_flux = dg.surface_integral.surface_flux # reshape face/normal arrays to have size = (num_points_on_face, num_faces_total). @@ -511,10 +509,10 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, # https://github.com/JuliaLang/julia/issues/36313#issuecomment-782336300. reshape_by_face(u) = Base.ReshapedArray(u, (num_pts_per_face, num_faces_total), ()) - u_face_values = reshape_by_face(u_face_values) - flux_face_values = reshape_by_face(flux_face_values) - Jf = reshape_by_face(Jf) - nxyzJ, xyzf = reshape_by_face.(nxyzJ), reshape_by_face.(xyzf) # broadcast over nxyzJ::NTuple{NDIMS,Matrix} + u_face_values = reshape_by_face(cache.u_face_values) + flux_face_values = reshape_by_face(cache.flux_face_values) + Jf = reshape_by_face(md.Jf) + nxyzJ, xyzf = reshape_by_face.(md.nxyzJ), reshape_by_face.(md.xyzf) # broadcast over nxyzJ::NTuple{NDIMS,Matrix} # loop through boundary faces, which correspond to columns of reshaped u_face_values, ... for f in mesh.boundary_faces[boundary_key] From d348ce32ba02ca82c69cc0711b62ea8c14dc14e1 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 16:21:17 -0500 Subject: [PATCH 18/35] replace some instances of `rd.Nfaces` with `StartUpDG.num_faces` --- src/solvers/dgmulti/dg.jl | 5 +++-- src/solvers/dgmulti/flux_differencing_gauss_sbp.jl | 10 ++++++---- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index f6c387d3db..f9c5578c18 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -459,8 +459,9 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, # reshape face/normal arrays to have size = (num_points_on_face, num_faces_total). # mesh.boundary_faces indexes into the columns of these face-reshaped arrays. - num_pts_per_face = rd.Nfq ÷ rd.Nfaces - num_faces_total = rd.Nfaces * md.num_elements + num_faces = StartUpDG.num_faces(rd.element_type) + num_pts_per_face = rd.Nfq ÷ num_faces + num_faces_total = num_faces * md.num_elements # This function was originally defined as # `reshape_by_face(u) = reshape(view(u, :), num_pts_per_face, num_faces_total)`. diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index b4542c527f..e2e18e1c00 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl @@ -61,7 +61,8 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, nnodes_1d = length(rq1D) # Permutation of indices in a tensor product form - indices = reshape(1:length(rd.rf), nnodes_1d, rd.Nfaces) + num_faces = StartUpDG.element_type(rd.element_type) + indices = reshape(1:length(rd.rf), nnodes_1d, num_faces) face_indices_tensor_product = zeros(Int, 2, nnodes_1d, ndims(rd.element_type)) for i in 1:nnodes_1d # loop over nodes in one face face_indices_tensor_product[:, i, 1] .= indices[i, 1:2] @@ -76,7 +77,7 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, return TensorProductGaussFaceOperator{2, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d, inv.(wq1D), rd.wf, face_indices_tensor_product, - nnodes_1d, rd.Nfaces) + nnodes_1d, num_faces) end # constructor for a 3D operator @@ -90,7 +91,8 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, nnodes_1d = length(rq1D) # Permutation of indices in a tensor product form - indices = reshape(1:length(rd.rf), nnodes_1d, nnodes_1d, rd.Nfaces) + num_faces = StartUpDG.num_faces(rd.element_type) + indices = reshape(1:length(rd.rf), nnodes_1d, nnodes_1d, num_faces) face_indices_tensor_product = zeros(Int, 2, nnodes_1d, nnodes_1d, ndims(rd.element_type)) for j in 1:nnodes_1d, i in 1:nnodes_1d # loop over nodes in one face face_indices_tensor_product[:, i, j, 1] .= indices[i, j, 1:2] @@ -106,7 +108,7 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, return TensorProductGaussFaceOperator{3, T_op, Tm, Tw, Tf, Ti}(interp_matrix_gauss_to_face_1d, inv.(wq1D), rd.wf, face_indices_tensor_product, - nnodes_1d, rd.Nfaces) + nnodes_1d, num_faces) end # specialize behavior of `mul_by!(A)` where `A isa TensorProductGaussFaceOperator)` From 4df3ed3cabd91edb7e3afb1b10080ddab5fcc4a2 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Mon, 8 May 2023 16:30:00 -0500 Subject: [PATCH 19/35] Update examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl Co-authored-by: Michael Schlottke-Lakemper --- examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index cd623109ed..64bff89ffc 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -82,7 +82,7 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() analysis_interval = 100 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg)) +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(solver)) alive_callback = AliveCallback(alive_interval=10) cfl = 0.5 From c0151b9f6eef790e494b6b7aef718c82382dbd3f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Mon, 8 May 2023 18:30:26 -0500 Subject: [PATCH 20/35] fix StartUpDG.num_faces call --- src/solvers/dgmulti/flux_differencing_gauss_sbp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index e2e18e1c00..c50da255d9 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl @@ -61,7 +61,7 @@ function TensorProductGaussFaceOperator(operator::AbstractGaussOperator, nnodes_1d = length(rq1D) # Permutation of indices in a tensor product form - num_faces = StartUpDG.element_type(rd.element_type) + num_faces = StartUpDG.num_faces(rd.element_type) indices = reshape(1:length(rd.rf), nnodes_1d, num_faces) face_indices_tensor_product = zeros(Int, 2, nnodes_1d, ndims(rd.element_type)) for i in 1:nnodes_1d # loop over nodes in one face From 16aaea7c539a2a6a0c0c3a3bf8dfefd60c0bceb4 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Tue, 9 May 2023 08:30:59 -0500 Subject: [PATCH 21/35] Update src/basic_types.jl Co-authored-by: Hendrik Ranocha --- src/basic_types.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/basic_types.jl b/src/basic_types.jl index b7bcede275..06f47518b0 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -84,7 +84,7 @@ end @inline function (::BoundaryConditionDoNothing)( u_inner, outward_direction::AbstractVector, x, t, surface_flux, equations) - # this should reduce to `flux(u_inner, outward_direction, equations)` for a consistent symmetric flux. + # this should reduce to `flux(u_inner, outward_direction, equations)` for a consistent flux. return surface_flux(u_inner, u_inner, outward_direction, equations) end From 2905886aa543d8b2ee44ae8c62544b867d7cb4ce Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Tue, 9 May 2023 08:32:13 -0500 Subject: [PATCH 22/35] Update examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl Co-authored-by: Hendrik Ranocha --- examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index 64bff89ffc..6afaa373c1 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -76,7 +76,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, .075) +tspan = (0.0, 0.075) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() From e9dfa39feaadf86c0081a338c886d2e8d9ebbe4a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 9 May 2023 08:35:23 -0500 Subject: [PATCH 23/35] update test elixir --- examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl | 6 +++--- test/test_dgmulti_2d.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index 64bff89ffc..18a183d006 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -1,6 +1,7 @@ using OrdinaryDiffEq using Trixi +using LinearAlgebra: norm, dot # for use in the MHD boundary condition ############################################################################### # semidiscretization of the compressible ideal GLM-MHD equations @@ -31,7 +32,7 @@ initial_condition = initial_condition_perturbation surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell) volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell) -solver = DGMulti(polydeg=3, element_type = Quad(), approximation_type = SBP(), +solver = DGMulti(polydeg=3, element_type = Quad(), approximation_type = GaussSBP(), surface_integral = SurfaceIntegralWeakForm(surface_flux), volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) @@ -46,7 +47,6 @@ mesh = DGMultiMesh(solver, cells_per_dimension; periodicity=(false, false), is_o # Create a "reflective-like" boundary condition by mirroring the velocity but leaving the magnetic field alone. # Note that this boundary condition is probably not entropy stable. -using LinearAlgebra: norm, dot function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector, x, t, surface_flux_function, equations::IdealGlmMhdEquations2D) @@ -90,7 +90,7 @@ stepsize_callback = StepsizeCallback(cfl=cfl) glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) callbacks = CallbackSet(summary_callback, - analysis_callback, + #analysis_callback, alive_callback, stepsize_callback, glm_speed_callback) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 9730486af4..6477a78cb3 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -294,7 +294,7 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "elixir_mhd_reflective_BCs.jl (Quad)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_reflective_BCs.jl"), - cells_per_dimension = 4, element_type = Quad(), approximation_type = GaussSBP(), + cells_per_dimension = 4, l2 = [0.0036019562526885497, 0.0017340971255526598, 0.008375221676923516, 0.0, 0.028596802653952254, 0.0018573697892179616, 0.002080779894047767, 0.0, 5.3012597624420295e-5], linf = [0.01692598382370969, 0.00936965952976002, 0.04145170727841342, 0.0, 0.11569901084179435, 0.009849648257862952, 0.011417088537134745, 0.0, 0.0002992621756978016] ) From edf028d9e6e9ba9e3bd550f70bca3d3270a67b0f Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Tue, 9 May 2023 08:36:22 -0500 Subject: [PATCH 24/35] Update src/solvers/dgmulti/dg.jl Co-authored-by: Hendrik Ranocha --- src/solvers/dgmulti/dg.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index f9c5578c18..aef8c080e9 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -428,8 +428,7 @@ end # do nothing for periodic (default) boundary conditions calc_boundary_flux!(cache, t, boundary_conditions::BoundaryConditionPeriodic, - have_nonconservative_terms, - mesh, equations, dg::DGMulti) = nothing + mesh, have_nonconservative_terms, equations, dg::DGMulti) = nothing # "lispy tuple programming" instead of for loop for type stability function calc_boundary_flux!(cache, t, boundary_conditions, have_nonconservative_terms, From 44689d593e16c87f3d409981ee43bdd146ab982d Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 9 May 2023 08:53:13 -0500 Subject: [PATCH 25/35] fix calc_boundary_flux! signature --- src/solvers/dgmulti/dg.jl | 25 +++++++++---------- src/solvers/dgmulti/flux_differencing.jl | 9 +++---- .../dgmulti/flux_differencing_gauss_sbp.jl | 5 ++-- 3 files changed, 18 insertions(+), 21 deletions(-) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index aef8c080e9..21f74f7fd8 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -431,24 +431,24 @@ calc_boundary_flux!(cache, t, boundary_conditions::BoundaryConditionPeriodic, mesh, have_nonconservative_terms, equations, dg::DGMulti) = nothing # "lispy tuple programming" instead of for loop for type stability -function calc_boundary_flux!(cache, t, boundary_conditions, have_nonconservative_terms, - mesh, equations, dg::DGMulti) +function calc_boundary_flux!(cache, t, boundary_conditions, mesh, + have_nonconservative_terms, equations, dg::DGMulti) + # peel off first boundary condition calc_single_boundary_flux!(cache, t, first(boundary_conditions), first(keys(boundary_conditions)), - have_nonconservative_terms, mesh, equations, dg) + mesh, have_nonconservative_terms, equations, dg) # recurse on the remainder of the boundary conditions calc_boundary_flux!(cache, t, Base.tail(boundary_conditions), - have_nonconservative_terms, mesh, equations, dg) + mesh, have_nonconservative_terms, equations, dg) end # terminate recursion calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple{(),Tuple{}}, - have_nonconservative_terms, mesh, equations, dg::DGMulti) = nothing + mesh, have_nonconservative_terms, equations, dg::DGMulti) = nothing -function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, - have_nonconservative_terms::False, - mesh, equations, dg::DGMulti{NDIMS}) where {NDIMS} +function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, mesh, + have_nonconservative_terms::False, equations, dg::DGMulti{NDIMS}) where {NDIMS} rd = dg.basis md = mesh.md @@ -489,9 +489,8 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, # However, we don't have to re-reshape, since cache.flux_face_values still retains its original shape. end -function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, - have_nonconservative_terms::True, - mesh, equations, dg::DGMulti{NDIMS}) where {NDIMS} +function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, mesh, + have_nonconservative_terms::True, equations, dg::DGMulti{NDIMS}) where {NDIMS} rd = dg.basis md = mesh.md @@ -619,8 +618,8 @@ function rhs!(du, u, t, mesh, equations, have_nonconservative_terms(equations), equations, dg) @trixi_timeit timer() "boundary flux" calc_boundary_flux!( - cache, t, boundary_conditions, - have_nonconservative_terms(equations), mesh, equations, dg) + cache, t, boundary_conditions, mesh, + have_nonconservative_terms(equations), equations, dg) @trixi_timeit timer() "surface integral" calc_surface_integral!( du, u, dg.surface_integral, mesh, equations, dg, cache) diff --git a/src/solvers/dgmulti/flux_differencing.jl b/src/solvers/dgmulti/flux_differencing.jl index b541941388..f511694c76 100644 --- a/src/solvers/dgmulti/flux_differencing.jl +++ b/src/solvers/dgmulti/flux_differencing.jl @@ -594,9 +594,8 @@ function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions: have_nonconservative_terms(equations), equations, dg) - @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, - have_nonconservative_terms(equations), - mesh, equations, dg) + @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, mesh, + have_nonconservative_terms(equations), equations, dg) @trixi_timeit timer() "surface integral" calc_surface_integral!(du, u, dg.surface_integral, mesh, equations, dg, cache) @@ -631,8 +630,8 @@ function rhs!(du, u, t, mesh, equations, have_nonconservative_terms(equations), equations, dg) @trixi_timeit timer() "boundary flux" calc_boundary_flux!( - cache, t, boundary_conditions, have_nonconservative_terms(equations), - mesh, equations, dg) + cache, t, boundary_conditions, mesh, + have_nonconservative_terms(equations), equations, dg) @trixi_timeit timer() "surface integral" calc_surface_integral!( du, u, dg.surface_integral, mesh, equations, dg, cache) diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl index c50da255d9..ca2666f218 100644 --- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl +++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl @@ -509,9 +509,8 @@ function rhs!(du, u, t, mesh, equations, initial_condition, boundary_conditions: have_nonconservative_terms(equations), equations, dg) - @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, - have_nonconservative_terms(equations), - mesh, equations, dg) + @trixi_timeit timer() "boundary flux" calc_boundary_flux!(cache, t, boundary_conditions, mesh, + have_nonconservative_terms(equations), equations, dg) # `du` is stored at Gauss nodes here @trixi_timeit timer() "surface integral" calc_surface_integral!(du, u, dg.surface_integral, From c0bc470f312e32c5a63cc73cc84aa44c9f665879 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 9 May 2023 10:25:55 -0500 Subject: [PATCH 26/35] switch to dispatch for Dirichlet/DoNothing BCs when using noncons flux --- .../dgmulti_2d/elixir_mhd_reflective_BCs.jl | 9 +++++++++ src/basic_types.jl | 14 ++++++++++++-- src/equations/equations.jl | 17 +++++++++++++++-- src/equations/ideal_glm_mhd_2d.jl | 5 ----- src/solvers/dgmulti/dg.jl | 5 +++-- 5 files changed, 39 insertions(+), 11 deletions(-) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index 637b413f46..3cadf95b0b 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -65,6 +65,15 @@ function boundary_condition_velocity_slip_wall(u_inner, normal_direction::Abstra return surface_flux_function(u_inner, u_mirror, normal, equations) * norm_ end +# we need to define two versions of boundary conditions when using non-conservative terms +@inline function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector, + normal_direction_avg::AbstractVector, x, t, + nonconservative_surface_flux, + equations::IdealGlmMhdEquations2D) + return boundary_condition_velocity_slip_wall(u_inner, normal_direction, x, t, + nonconservative_surface_flux, equations) +end + boundary_conditions = (; x_neg=boundary_condition_velocity_slip_wall, x_pos=boundary_condition_velocity_slip_wall, y_neg=boundary_condition_do_nothing, diff --git a/src/basic_types.jl b/src/basic_types.jl index 06f47518b0..c9d903fc3c 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -76,18 +76,28 @@ struct BoundaryConditionDoNothing end # This version can be called by hyperbolic solvers on logically Cartesian meshes @inline function (::BoundaryConditionDoNothing)( u_inner, orientation_or_normal_direction, direction::Integer, x, t, surface_flux, equations) + return flux(u_inner, orientation_or_normal_direction, equations) # TODO: should this be switched to `surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)` for consistency? end # This version can be called by hyperbolic solvers on unstructured, curved meshes -@inline function (::BoundaryConditionDoNothing)( - u_inner, outward_direction::AbstractVector, x, t, surface_flux, equations) +@inline function (::BoundaryConditionDoNothing)(u_inner, outward_direction::AbstractVector, + x, t, surface_flux, equations) # this should reduce to `flux(u_inner, outward_direction, equations)` for a consistent flux. return surface_flux(u_inner, u_inner, outward_direction, equations) end +# This version is called when passing in a non-conservative surface flux +@inline function (::BoundaryConditionDoNothing)(u_inner, outward_direction::AbstractVector, + outward_direction_avg::AbstractVector, + x, t, nonconservative_surface_flux, equations) + + return nonconservative_surface_flux(u_inner, u_inner, outward_direction, + outward_direction_avg, equations) +end + # This version can be called by parabolic solvers @inline function (::BoundaryConditionDoNothing)(inner_flux_or_state, other_args...) return inner_flux_or_state diff --git a/src/equations/equations.jl b/src/equations/equations.jl index c1669531de..16aba6dd66 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -181,6 +181,19 @@ end return flux end +@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, + normal_direction::AbstractVector, + normal_direction_avg::AbstractVector, + x, t, nonconservative_surface_flux, equations) + # get the external value of the solution + u_boundary = boundary_condition.boundary_value_function(x, t, equations) + + # Calculate boundary flux + flux = nonconservative_surface_flux(u_inner, u_boundary, normal_direction, normal_direction_avg, equations) + + return flux +end + # operator types used for dispatch on parabolic boundary fluxes struct Gradient end struct Divergence end @@ -328,7 +341,7 @@ include("compressible_euler_multicomponent_2d.jl") eachcomponent(equations::AbstractCompressibleEulerMulticomponentEquations) Return an iterator over the indices that specify the location in relevant data structures -for the components in `AbstractCompressibleEulerMulticomponentEquations`. +for the components in `AbstractCompressibleEulerMulticomponentEquations`. In particular, not the components themselves are returned. """ @inline eachcomponent(equations::AbstractCompressibleEulerMulticomponentEquations) = Base.OneTo(ncomponents(equations)) @@ -350,7 +363,7 @@ include("ideal_glm_mhd_multicomponent_2d.jl") eachcomponent(equations::AbstractIdealGlmMhdMulticomponentEquations) Return an iterator over the indices that specify the location in relevant data structures -for the components in `AbstractIdealGlmMhdMulticomponentEquations`. +for the components in `AbstractIdealGlmMhdMulticomponentEquations`. In particular, not the components themselves are returned. """ @inline eachcomponent(equations::AbstractIdealGlmMhdMulticomponentEquations) = Base.OneTo(ncomponents(equations)) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index f51e8f652d..c19273737e 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -238,11 +238,6 @@ terms. return f end -# TODO: is there a way to do this without manually dispatching for every non-conservative flux? We use this for the evaluation of boundary conditions. -# if only one normal is specified, use that for `normal_direction_average` -@inline flux_nonconservative_powell(u_ll, u_rr, normal_direction::AbstractVector, equations::IdealGlmMhdEquations2D) = - flux_nonconservative_powell(u_ll, u_rr, normal_direction, normal_direction, equations) - @inline function flux_nonconservative_powell(u_ll, u_rr, normal_direction_ll::AbstractVector, normal_direction_average::AbstractVector, diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 21f74f7fd8..46f37744b8 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -523,8 +523,9 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, cons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, surface_flux, equations) - noncons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, - nonconservative_flux, equations) + # the non-conservative flux involves two face normals + noncons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_normal, + face_coordinates, t, nonconservative_flux, equations) flux_face_values[i,f] = (cons_flux_at_face_node + 0.5 * noncons_flux_at_face_node) * Jf[i,f] From 5773509fcab5b0eca852239a6d5653f39d4e9a2a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 9 May 2023 11:18:07 -0500 Subject: [PATCH 27/35] fix nonconservative BC --- .../dgmulti_2d/elixir_mhd_reflective_BCs.jl | 20 +++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index 3cadf95b0b..a22038757c 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -70,8 +70,24 @@ end normal_direction_avg::AbstractVector, x, t, nonconservative_surface_flux, equations::IdealGlmMhdEquations2D) - return boundary_condition_velocity_slip_wall(u_inner, normal_direction, x, t, - nonconservative_surface_flux, equations) + # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later + norm_ = norm(normal_direction) + normal = normal_direction / norm_ + norm_avg_ = norm(normal_direction_avg) + normal_avg = normal_direction_avg / norm_avg_ # assume both have the same norm + + # average the magnitudes of the two normals in case they differ + norm_ = 0.5 * (norm_ + norm_avg_) + + # compute the primitive variables + rho, v1, v2, v3, p, B1, B2, B3, psi = cons2prim(u_inner, equations) + + v_normal = dot(normal, SVector(v1, v2)) + u_mirror = prim2cons(SVector(rho, v1 - 2 * v_normal * normal[1], + v2 - 2 * v_normal * normal[2], + v3, p, B1, B2, B3, psi), equations) + + return nonconservative_surface_flux(u_inner, u_mirror, normal, normal_avg, equations) * norm_ end boundary_conditions = (; x_neg=boundary_condition_velocity_slip_wall, From df9a71abe1e87c50a457b6ca6af4bdce63f6792a Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 9 May 2023 11:53:06 -0500 Subject: [PATCH 28/35] fix type ambiguity --- src/equations/equations.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 16aba6dd66..b07557c83d 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -171,7 +171,8 @@ end @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal_direction::AbstractVector, x, t, - surface_flux_function, equations) + surface_flux_function, + equations::AbstractEquations) # get the external value of the solution u_boundary = boundary_condition.boundary_value_function(x, t, equations) @@ -184,7 +185,8 @@ end @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal_direction::AbstractVector, normal_direction_avg::AbstractVector, - x, t, nonconservative_surface_flux, equations) + x, t, nonconservative_surface_flux, + equations::AbstractEquations) # get the external value of the solution u_boundary = boundary_condition.boundary_value_function(x, t, equations) From 6a184afb3fa4663e4fcd1e37a8780d7cdcd6121f Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Tue, 9 May 2023 13:33:19 -0500 Subject: [PATCH 29/35] fix type ambiguity by redesigning nonconservative BC signature --- examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl | 11 +++-------- src/basic_types.jl | 6 +++--- src/equations/equations.jl | 10 +++++----- src/solvers/dgmulti/dg.jl | 7 ++++--- 4 files changed, 15 insertions(+), 19 deletions(-) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index a22038757c..d6f764eeed 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -65,19 +65,14 @@ function boundary_condition_velocity_slip_wall(u_inner, normal_direction::Abstra return surface_flux_function(u_inner, u_mirror, normal, equations) * norm_ end -# we need to define two versions of boundary conditions when using non-conservative terms +# we need to define two types of boundary conditions when using non-conservative terms @inline function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector, - normal_direction_avg::AbstractVector, x, t, + x, t, flux_is_nonconservative::Trixi.True, nonconservative_surface_flux, equations::IdealGlmMhdEquations2D) # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later norm_ = norm(normal_direction) normal = normal_direction / norm_ - norm_avg_ = norm(normal_direction_avg) - normal_avg = normal_direction_avg / norm_avg_ # assume both have the same norm - - # average the magnitudes of the two normals in case they differ - norm_ = 0.5 * (norm_ + norm_avg_) # compute the primitive variables rho, v1, v2, v3, p, B1, B2, B3, psi = cons2prim(u_inner, equations) @@ -87,7 +82,7 @@ end v2 - 2 * v_normal * normal[2], v3, p, B1, B2, B3, psi), equations) - return nonconservative_surface_flux(u_inner, u_mirror, normal, normal_avg, equations) * norm_ + return nonconservative_surface_flux(u_inner, u_mirror, normal, normal, equations) * norm_ end boundary_conditions = (; x_neg=boundary_condition_velocity_slip_wall, diff --git a/src/basic_types.jl b/src/basic_types.jl index c9d903fc3c..04e0bf995e 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -91,11 +91,11 @@ end # This version is called when passing in a non-conservative surface flux @inline function (::BoundaryConditionDoNothing)(u_inner, outward_direction::AbstractVector, - outward_direction_avg::AbstractVector, - x, t, nonconservative_surface_flux, equations) + x, t, using_nonconservative_flux::True, + nonconservative_surface_flux, equations) return nonconservative_surface_flux(u_inner, u_inner, outward_direction, - outward_direction_avg, equations) + outward_direction, equations) end # This version can be called by parabolic solvers diff --git a/src/equations/equations.jl b/src/equations/equations.jl index b07557c83d..afc7b61fcb 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -172,7 +172,7 @@ end normal_direction::AbstractVector, x, t, surface_flux_function, - equations::AbstractEquations) + equations) # get the external value of the solution u_boundary = boundary_condition.boundary_value_function(x, t, equations) @@ -184,14 +184,14 @@ end @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal_direction::AbstractVector, - normal_direction_avg::AbstractVector, - x, t, nonconservative_surface_flux, - equations::AbstractEquations) + x, t, using_nonconservative_flux::True, + nonconservative_surface_flux, + equations) # get the external value of the solution u_boundary = boundary_condition.boundary_value_function(x, t, equations) # Calculate boundary flux - flux = nonconservative_surface_flux(u_inner, u_boundary, normal_direction, normal_direction_avg, equations) + flux = nonconservative_surface_flux(u_inner, u_boundary, normal_direction, normal_direction, equations) return flux end diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 46f37744b8..2eee55e14f 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -523,9 +523,10 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, cons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, surface_flux, equations) - # the non-conservative flux involves two face normals - noncons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_normal, - face_coordinates, t, nonconservative_flux, equations) + # specify that we are using a non-conservative flux here + using_nonconservative_flux = True() + noncons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, + using_nonconservative_flux, nonconservative_flux, equations) flux_face_values[i,f] = (cons_flux_at_face_node + 0.5 * noncons_flux_at_face_node) * Jf[i,f] From 7ea24901ad6b725e9a7947346e205e69aba24352 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Wed, 10 May 2023 08:15:51 -0500 Subject: [PATCH 30/35] Update src/basic_types.jl Co-authored-by: Hendrik Ranocha --- src/basic_types.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/basic_types.jl b/src/basic_types.jl index 04e0bf995e..681f785212 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -85,8 +85,7 @@ end @inline function (::BoundaryConditionDoNothing)(u_inner, outward_direction::AbstractVector, x, t, surface_flux, equations) - # this should reduce to `flux(u_inner, outward_direction, equations)` for a consistent flux. - return surface_flux(u_inner, u_inner, outward_direction, equations) + return flux(u_inner, outward_direction, equations) end # This version is called when passing in a non-conservative surface flux From 15e1a5e6a81aaeae48af852b10080de0614da04f Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Wed, 10 May 2023 08:16:21 -0500 Subject: [PATCH 31/35] Update examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl Co-authored-by: Hendrik Ranocha --- examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index d6f764eeed..654bb274f3 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -110,7 +110,7 @@ stepsize_callback = StepsizeCallback(cfl=cfl) glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl) callbacks = CallbackSet(summary_callback, - #analysis_callback, + analysis_callback, alive_callback, stepsize_callback, glm_speed_callback) From 24d355b76385fa1f3b15e9be67cf5a89c2ddfbf2 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Wed, 10 May 2023 08:59:35 -0500 Subject: [PATCH 32/35] Update src/basic_types.jl Co-authored-by: Hendrik Ranocha --- src/basic_types.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/basic_types.jl b/src/basic_types.jl index 681f785212..0049656354 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -78,7 +78,6 @@ struct BoundaryConditionDoNothing end u_inner, orientation_or_normal_direction, direction::Integer, x, t, surface_flux, equations) return flux(u_inner, orientation_or_normal_direction, equations) - # TODO: should this be switched to `surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)` for consistency? end # This version can be called by hyperbolic solvers on unstructured, curved meshes From faae89716310d4b41e4a2abf1d49976d536a5a23 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Wed, 10 May 2023 09:29:21 -0500 Subject: [PATCH 33/35] make nonconservative BCs consistent with rest of Trixi --- .../dgmulti_2d/elixir_mhd_reflective_BCs.jl | 20 ------------------- src/solvers/dgmulti/dg.jl | 14 ++++++++----- test/test_dgmulti_2d.jl | 4 ++-- 3 files changed, 11 insertions(+), 27 deletions(-) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl index 654bb274f3..a1351cf824 100644 --- a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl +++ b/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl @@ -65,26 +65,6 @@ function boundary_condition_velocity_slip_wall(u_inner, normal_direction::Abstra return surface_flux_function(u_inner, u_mirror, normal, equations) * norm_ end -# we need to define two types of boundary conditions when using non-conservative terms -@inline function boundary_condition_velocity_slip_wall(u_inner, normal_direction::AbstractVector, - x, t, flux_is_nonconservative::Trixi.True, - nonconservative_surface_flux, - equations::IdealGlmMhdEquations2D) - # Normalize the vector without using `normalize` since we need to multiply by the `norm_` later - norm_ = norm(normal_direction) - normal = normal_direction / norm_ - - # compute the primitive variables - rho, v1, v2, v3, p, B1, B2, B3, psi = cons2prim(u_inner, equations) - - v_normal = dot(normal, SVector(v1, v2)) - u_mirror = prim2cons(SVector(rho, v1 - 2 * v_normal * normal[1], - v2 - 2 * v_normal * normal[2], - v3, p, B1, B2, B3, psi), equations) - - return nonconservative_surface_flux(u_inner, u_mirror, normal, normal, equations) * norm_ -end - boundary_conditions = (; x_neg=boundary_condition_velocity_slip_wall, x_pos=boundary_condition_velocity_slip_wall, y_neg=boundary_condition_do_nothing, diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 2eee55e14f..5d087f0deb 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -519,14 +519,18 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, face_normal = SVector{NDIMS}(getindex.(nxyzJ, i, f)) / Jf[i,f] face_coordinates = SVector{NDIMS}(getindex.(xyzf, i, f)) - # compute conservative and non-conservative parts separately + # Compute conservative and non-conservative fluxes separately. + # This imposes boundary conditions on the conservative part of the flux. cons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, surface_flux, equations) - # specify that we are using a non-conservative flux here - using_nonconservative_flux = True() - noncons_flux_at_face_node = boundary_condition(u_face_values[i,f], face_normal, face_coordinates, t, - using_nonconservative_flux, nonconservative_flux, equations) + # Compute pointwise nonconservative numerical flux at the boundary. + # In general, nonconservative fluxes can depend on both the contravariant + # vectors (normal direction) at the current node and the averaged ones. + # However, there is only one `face_normal` at boundaries, which we pass in twice. + # Note: This does not set any type of boundary condition for the nonconservative term + noncons_flux_at_face_node = nonconservative_flux(u_face_values[i,f], u_face_values[i,f], + face_normal, face_normal, equations) flux_face_values[i,f] = (cons_flux_at_face_node + 0.5 * noncons_flux_at_face_node) * Jf[i,f] diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 6477a78cb3..3786f879fc 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -295,8 +295,8 @@ isdir(outdir) && rm(outdir, recursive=true) @trixi_testset "elixir_mhd_reflective_BCs.jl (Quad)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_reflective_BCs.jl"), cells_per_dimension = 4, - l2 = [0.0036019562526885497, 0.0017340971255526598, 0.008375221676923516, 0.0, 0.028596802653952254, 0.0018573697892179616, 0.002080779894047767, 0.0, 5.3012597624420295e-5], - linf = [0.01692598382370969, 0.00936965952976002, 0.04145170727841342, 0.0, 0.11569901084179435, 0.009849648257862952, 0.011417088537134745, 0.0, 0.0002992621756978016] + l2 = [0.0036019536614619687, 0.001734097206958611, 0.008375221008997178, 0.0, 0.028596796602124414, 0.0018573693138866614, 0.0020807798141551166, 0.0, 5.301188920230166e-5], + linf = [0.01692601228199253, 0.009369662298436778, 0.04145169295835428, 0.0, 0.11569908670112738, 0.00984964453299233, 0.01141708032148614, 0.0, 0.0002992631411931389] ) end From fd8d6619c58386ce9a87a4b3e3c46a714f563fd7 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 12 May 2023 10:46:36 -0500 Subject: [PATCH 34/35] renaming --- ...ir_mhd_reflective_BCs.jl => elixir_mhd_reflective_wall.jl} | 0 test/test_dgmulti_2d.jl | 4 ++-- 2 files changed, 2 insertions(+), 2 deletions(-) rename examples/dgmulti_2d/{elixir_mhd_reflective_BCs.jl => elixir_mhd_reflective_wall.jl} (100%) diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl b/examples/dgmulti_2d/elixir_mhd_reflective_wall.jl similarity index 100% rename from examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl rename to examples/dgmulti_2d/elixir_mhd_reflective_wall.jl diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 3786f879fc..0c10a17642 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -292,8 +292,8 @@ isdir(outdir) && rm(outdir, recursive=true) ) end - @trixi_testset "elixir_mhd_reflective_BCs.jl (Quad)" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_reflective_BCs.jl"), + @trixi_testset "elixir_mhd_reflective_wall.jl (Quad)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_reflective_wall.jl"), cells_per_dimension = 4, l2 = [0.0036019536614619687, 0.001734097206958611, 0.008375221008997178, 0.0, 0.028596796602124414, 0.0018573693138866614, 0.0020807798141551166, 0.0, 5.301188920230166e-5], linf = [0.01692601228199253, 0.009369662298436778, 0.04145169295835428, 0.0, 0.11569908670112738, 0.00984964453299233, 0.01141708032148614, 0.0, 0.0002992631411931389] From 091c5b5374d8cd3345ce61ebbb53745182e96ac0 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 12 May 2023 10:58:52 -0500 Subject: [PATCH 35/35] deleting unused boundary condition implementations --- src/basic_types.jl | 9 --------- src/equations/equations.jl | 14 -------------- 2 files changed, 23 deletions(-) diff --git a/src/basic_types.jl b/src/basic_types.jl index 0049656354..4539e26dea 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -87,15 +87,6 @@ end return flux(u_inner, outward_direction, equations) end -# This version is called when passing in a non-conservative surface flux -@inline function (::BoundaryConditionDoNothing)(u_inner, outward_direction::AbstractVector, - x, t, using_nonconservative_flux::True, - nonconservative_surface_flux, equations) - - return nonconservative_surface_flux(u_inner, u_inner, outward_direction, - outward_direction, equations) -end - # This version can be called by parabolic solvers @inline function (::BoundaryConditionDoNothing)(inner_flux_or_state, other_args...) return inner_flux_or_state diff --git a/src/equations/equations.jl b/src/equations/equations.jl index afc7b61fcb..e44270737e 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -182,20 +182,6 @@ end return flux end -@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, - normal_direction::AbstractVector, - x, t, using_nonconservative_flux::True, - nonconservative_surface_flux, - equations) - # get the external value of the solution - u_boundary = boundary_condition.boundary_value_function(x, t, equations) - - # Calculate boundary flux - flux = nonconservative_surface_flux(u_inner, u_boundary, normal_direction, normal_direction, equations) - - return flux -end - # operator types used for dispatch on parabolic boundary fluxes struct Gradient end struct Divergence end