Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

specialize calc_boundary_flux! for nonconservative terms for DGMulti #1431

Merged
merged 42 commits into from
May 14, 2023
Merged
Show file tree
Hide file tree
Changes from 16 commits
Commits
Show all changes
42 commits
Select commit Hold shift + click to select a range
ba70d67
minor change for consistency
jlchan May 8, 2023
51c2a5e
formatting
jlchan May 8, 2023
059c0e8
add nonconservative terms to DGMulti `calc_boundary_flux!`
jlchan May 8, 2023
7d820b9
add noncon boundary flux
jlchan May 8, 2023
e2d1432
fix dropped dg.surface_flux
jlchan May 8, 2023
3fc2254
Merge branch 'main' into jc/DGMulti_noncon_BCs
jlchan May 8, 2023
f415e38
Merge remote-tracking branch 'Trixi_fork/jc/DGMulti_noncon_BCs' into …
jlchan May 8, 2023
cafb3f1
formatting
jlchan May 8, 2023
8987e9a
clean up noncons BCs
jlchan May 8, 2023
92ac9be
adding specialization of nonconservative Powell flux for BCs
jlchan May 8, 2023
ba9801d
fix BoundaryConditionDoNothing for nonconservative terms
jlchan May 8, 2023
22899b4
add elixir
jlchan May 8, 2023
24a4f1e
add test
jlchan May 8, 2023
0596c32
comment
jlchan May 8, 2023
e05b731
importing norm
jlchan May 8, 2023
f7131e4
import dot as well
jlchan May 8, 2023
aefa0cd
adding forgotten analysis callback
jlchan May 8, 2023
271399a
Update src/solvers/dgmulti/dg.jl
jlchan May 8, 2023
fef30d1
Merge remote-tracking branch 'Trixi_fork/jc/DGMulti_noncon_BCs' into …
jlchan May 8, 2023
383eb24
remove some name-based type instabilities
jlchan May 8, 2023
d348ce3
replace some instances of `rd.Nfaces` with `StartUpDG.num_faces`
jlchan May 8, 2023
4df3ed3
Update examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl
jlchan May 8, 2023
c0151b9
fix StartUpDG.num_faces call
jlchan May 8, 2023
16aaea7
Update src/basic_types.jl
jlchan May 9, 2023
d3dd0b0
Merge branch 'main' into jc/DGMulti_noncon_BCs
jlchan May 9, 2023
2905886
Update examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl
jlchan May 9, 2023
e9dfa39
update test elixir
jlchan May 9, 2023
c6e7c95
Merge remote-tracking branch 'Trixi_fork/jc/DGMulti_noncon_BCs' into …
jlchan May 9, 2023
edf028d
Update src/solvers/dgmulti/dg.jl
jlchan May 9, 2023
44689d5
fix calc_boundary_flux! signature
jlchan May 9, 2023
c0bc470
switch to dispatch for Dirichlet/DoNothing BCs when using noncons flux
jlchan May 9, 2023
5773509
fix nonconservative BC
jlchan May 9, 2023
df9a71a
fix type ambiguity
jlchan May 9, 2023
6a184af
fix type ambiguity by redesigning nonconservative BC signature
jlchan May 9, 2023
7ea2490
Update src/basic_types.jl
jlchan May 10, 2023
15e1a5e
Update examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl
jlchan May 10, 2023
24d355b
Update src/basic_types.jl
jlchan May 10, 2023
faae897
make nonconservative BCs consistent with rest of Trixi
jlchan May 10, 2023
c9c272e
Merge branch 'main' into jc/DGMulti_noncon_BCs
jlchan May 11, 2023
fd8d661
renaming
jlchan May 12, 2023
30d040b
Merge branch 'main' into jc/DGMulti_noncon_BCs
jlchan May 12, 2023
091c5b5
deleting unused boundary condition implementations
jlchan May 12, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
102 changes: 102 additions & 0 deletions examples/dgmulti_2d/elixir_mhd_reflective_BCs.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@

jlchan marked this conversation as resolved.
Show resolved Hide resolved
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)

# 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
jlchan marked this conversation as resolved.
Show resolved Hide resolved
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)
jlchan marked this conversation as resolved.
Show resolved Hide resolved
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
7 changes: 5 additions & 2 deletions src/basic_types.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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?
jlchan marked this conversation as resolved.
Show resolved Hide resolved
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.
jlchan marked this conversation as resolved.
Show resolved Hide resolved
return surface_flux(u_inner, u_inner, outward_direction, equations)
jlchan marked this conversation as resolved.
Show resolved Hide resolved
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()

Expand Down
2 changes: 1 addition & 1 deletion src/callbacks_step/glm_speed_dg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
5 changes: 5 additions & 0 deletions src/equations/ideal_glm_mhd_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
jlchan marked this conversation as resolved.
Show resolved Hide resolved

@inline function flux_nonconservative_powell(u_ll, u_rr,
normal_direction_ll::AbstractVector,
normal_direction_average::AbstractVector,
Expand Down
64 changes: 59 additions & 5 deletions src/solvers/dgmulti/dg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
jlchan marked this conversation as resolved.
Show resolved Hide resolved

# "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
Expand Down Expand Up @@ -485,6 +489,55 @@ 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 = 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 ÷ rd.Nfaces
num_faces_total = rd.Nfaces * md.num_elements
jlchan marked this conversation as resolved.
Show resolved Hide resolved

# 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}
jlchan marked this conversation as resolved.
Show resolved Hide resolved

# 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
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)

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)
Expand Down Expand Up @@ -568,7 +621,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)
Expand Down
4 changes: 3 additions & 1 deletion src/solvers/dgmulti/flux_differencing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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)
Expand Down
1 change: 1 addition & 0 deletions src/solvers/dgmulti/flux_differencing_gauss_sbp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/solvers/dgsem_tree/dg_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 8 additions & 0 deletions test/test_dgmulti_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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(),
jlchan marked this conversation as resolved.
Show resolved Hide resolved
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(),
Expand Down