From 8d8619c975380296fb03322811215142120a04ee Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Sun, 14 May 2023 00:28:50 -0500 Subject: [PATCH] specialize `calc_boundary_flux!` for nonconservative terms for `DGMulti` (#1431) * minor change for consistency * formatting * add nonconservative terms to DGMulti `calc_boundary_flux!` * add noncon boundary flux * fix dropped dg.surface_flux * formatting * clean up noncons BCs * adding specialization of nonconservative Powell flux for BCs * fix BoundaryConditionDoNothing for nonconservative terms * add elixir * add test * comment * importing norm * import dot as well * adding forgotten analysis callback * Update src/solvers/dgmulti/dg.jl Co-authored-by: Michael Schlottke-Lakemper * remove some name-based type instabilities * replace some instances of `rd.Nfaces` with `StartUpDG.num_faces` * Update examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl Co-authored-by: Michael Schlottke-Lakemper * fix StartUpDG.num_faces call * Update src/basic_types.jl Co-authored-by: Hendrik Ranocha * Update examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl Co-authored-by: Hendrik Ranocha * update test elixir * Update src/solvers/dgmulti/dg.jl Co-authored-by: Hendrik Ranocha * fix calc_boundary_flux! signature * switch to dispatch for Dirichlet/DoNothing BCs when using noncons flux * fix nonconservative BC * fix type ambiguity * fix type ambiguity by redesigning nonconservative BC signature * Update src/basic_types.jl Co-authored-by: Hendrik Ranocha * Update examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl Co-authored-by: Hendrik Ranocha * Update src/basic_types.jl Co-authored-by: Hendrik Ranocha * make nonconservative BCs consistent with rest of Trixi * renaming * deleting unused boundary condition implementations --------- Co-authored-by: Michael Schlottke-Lakemper Co-authored-by: Hendrik Ranocha --- .../dgmulti_2d/elixir_mhd_reflective_wall.jl | 105 ++++++++++++++++++ src/basic_types.jl | 8 +- src/callbacks_step/glm_speed_dg.jl | 2 +- src/equations/equations.jl | 7 +- src/solvers/dgmulti/dg.jl | 77 +++++++++++-- src/solvers/dgmulti/flux_differencing.jl | 7 +- .../dgmulti/flux_differencing_gauss_sbp.jl | 14 ++- src/solvers/dgsem_tree/dg_2d.jl | 2 +- test/test_dgmulti_2d.jl | 8 ++ 9 files changed, 203 insertions(+), 27 deletions(-) create mode 100644 examples/dgmulti_2d/elixir_mhd_reflective_wall.jl diff --git a/examples/dgmulti_2d/elixir_mhd_reflective_wall.jl b/examples/dgmulti_2d/elixir_mhd_reflective_wall.jl new file mode 100644 index 0000000000..a1351cf824 --- /dev/null +++ b/examples/dgmulti_2d/elixir_mhd_reflective_wall.jl @@ -0,0 +1,105 @@ + +using OrdinaryDiffEq +using Trixi +using LinearAlgebra: norm, dot # for use in the MHD boundary condition + +############################################################################### +# 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 = GaussSBP(), + 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) + +# 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) + + # 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, 0.075) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(solver)) +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, + analysis_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 diff --git a/src/basic_types.jl b/src/basic_types.jl index ea7d542c5d..4539e26dea 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -76,12 +76,14 @@ 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) 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) + return flux(u_inner, outward_direction, equations) end @@ -89,7 +91,7 @@ end @inline function (::BoundaryConditionDoNothing)(inner_flux_or_state, other_args...) return inner_flux_or_state end - + """ boundary_condition_do_nothing = Trixi.BoundaryConditionDoNothing() 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 diff --git a/src/equations/equations.jl b/src/equations/equations.jl index c1669531de..e44270737e 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) # get the external value of the solution u_boundary = boundary_condition.boundary_value_function(x, t, equations) @@ -328,7 +329,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 +351,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/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index f9e30f8f87..5d087f0deb 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/src/solvers/dgmulti/dg.jl @@ -428,24 +428,27 @@ end # do nothing for periodic (default) boundary conditions calc_boundary_flux!(cache, t, boundary_conditions::BoundaryConditionPeriodic, - 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, 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)), - 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), mesh, equations, dg) + calc_boundary_flux!(cache, t, Base.tail(boundary_conditions), + mesh, have_nonconservative_terms, equations, dg) end # terminate recursion calc_boundary_flux!(cache, t, boundary_conditions::NamedTuple{(),Tuple{}}, - mesh, equations, dg::DGMulti) = nothing + mesh, have_nonconservative_terms, equations, dg::DGMulti) = nothing -function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key, - 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 @@ -455,8 +458,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)`. @@ -485,6 +489,58 @@ 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, mesh, + have_nonconservative_terms::True, equations, dg::DGMulti{NDIMS}) where {NDIMS} + + rd = dg.basis + md = mesh.md + 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. + 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)`. + # 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(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] + 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 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) + + # 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] + + 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) @@ -568,7 +624,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, 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 97905f1d0b..f511694c76 100644 --- a/src/solvers/dgmulti/flux_differencing.jl +++ b/src/solvers/dgmulti/flux_differencing.jl @@ -594,8 +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, - 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) @@ -630,7 +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, 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 dfafe4ff98..ca2666f218 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.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 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)` @@ -507,8 +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, - 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, 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 diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index d2eaeea57d..0c10a17642 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_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] + ) + 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(),