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.

Sign up for GitHub

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

WIP: Add flexibility to nonconservative BCs #2200

Draft
wants to merge 7 commits into
base: main
Choose a base branch
from
Draft

WIP: Add flexibility to nonconservative BCs #2200

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
e2e2c20
adding flexibility to noncons BCs
MarcoArtiano Dec 10, 2024
1e925fc
minor fix
MarcoArtiano Dec 10, 2024
31f1e0a
Dirichlet BC for noncons
MarcoArtiano Dec 10, 2024
eb1ccdd
formatting
MarcoArtiano Dec 10, 2024
e4bc181
Noncons Dirichlet BC
MarcoArtiano Dec 10, 2024
8c5d2fc
test_type.jl fix
MarcoArtiano Dec 10, 2024
c2b44ae
dg multi fix
MarcoArtiano Dec 10, 2024
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
9 changes: 5 additions & 4 deletions examples/dgmulti_2d/elixir_mhd_reflective_wall.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,9 +49,9 @@ mesh = DGMultiMesh(solver, cells_per_dimension; periodicity = (false, false),
# 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,
x, t, surface_flux_functions,
equations::IdealGlmMhdEquations2D)

surface_flux_function, nonconservative_flux_function = surface_flux_functions
# Normalize the vector without using `normalize` since we need to multiply by the `norm_` later
norm_ = norm(normal_direction)
normal = normal_direction / norm_
Expand All @@ -63,8 +63,9 @@ function boundary_condition_velocity_slip_wall(u_inner, normal_direction::Abstra
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_
flux = surface_flux_function(u_inner, u_mirror, normal, equations) * norm_
noncons = nonconservative_flux_function(u_inner, u_mirror, normal, equations) * norm_
return flux + 0.5f0 * noncons
end

boundary_conditions = (; x_neg = boundary_condition_velocity_slip_wall,
Expand Down
24 changes: 24 additions & 0 deletions src/basic_types.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,13 +74,37 @@ struct BoundaryConditionDoNothing end
return surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)
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_functions::Tuple,
equations)
surface_flux_function, nonconservative_flux_function = surface_flux_functions
return surface_flux_function(u_inner, u_inner, orientation_or_normal_direction,
equations) +
0.5f0 * nonconservative_flux_function(u_inner, 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)
return surface_flux(u_inner, u_inner, outward_direction, equations)
end

@inline function (::BoundaryConditionDoNothing)(u_inner,
outward_direction::AbstractVector,
x, t, surface_flux_functions::Tuple,
equations)
surface_flux_function, nonconservative_flux_function = surface_flux_functions

return surface_flux_function(u_inner, u_inner, outward_direction, equations) +
0.5f0 *
nonconservative_flux_function(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
Expand Down
48 changes: 48 additions & 0 deletions src/equations/equations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,6 +184,35 @@ end
return flux
end

# Dirichlet-type boundary condition for use with TreeMesh or StructuredMesh
@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
orientation_or_normal,
direction,
x, t,
surface_flux_functions::Tuple,
equations)
surface_flux_function, nonconservative_flux_function = surface_flux_functions

u_boundary = boundary_condition.boundary_value_function(x, t, equations)

# Calculate boundary flux
if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal,
equations)
noncons = nonconservative_flux_function(u_inner, u_boundary,
orientation_or_normal,
equations)
else # u_boundary is "left" of boundary, u_inner is "right" of boundary
flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal,
equations)
noncons = nonconservative_flux_function(u_boundary, u_inner,
orientation_or_normal,
equations)
end

return flux + 0.5f0 * noncons
end

# Dirichlet-type boundary condition for use with UnstructuredMesh2D
# Note: For unstructured we lose the concept of an "absolute direction"
@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
Expand All @@ -200,6 +229,25 @@ end
return flux
end

# Dirichlet-type boundary condition for equations with non-conservative terms for use with UnstructuredMesh2D
# Note: For unstructured we lose the concept of an "absolute direction"
@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner,
normal_direction::AbstractVector,
x, t,
surface_flux_functions::Tuple,
equations)
surface_flux_function, nonconservative_flux_function = surface_flux_functions

# get the external value of the solution
u_boundary = boundary_condition.boundary_value_function(x, t, equations)

# Calculate boundary flux
flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
noncons = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
equations)
return flux + 0.5f0 * noncons
end

# operator types used for dispatch on parabolic boundary fluxes
struct Gradient end
struct Divergence end
Expand Down
12 changes: 10 additions & 2 deletions src/equations/shallow_water_1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,9 +160,11 @@ For details see Section 9.2.5 of the book:
"""
@inline function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction,
x, t,
surface_flux_function,
surface_flux_functions,
equations::ShallowWaterEquations1D)

# The boundary conditions for the non-conservative term are identically 0 here.
surface_flux_function, nonconservative_flux_function = surface_flux_functions
# create the "external" boundary solution state
u_boundary = SVector(u_inner[1],
-u_inner[2],
Expand All @@ -172,12 +174,18 @@ For details see Section 9.2.5 of the book:
if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal,
equations)
noncons = nonconservative_flux_function(u_inner, u_boundary,
orientation_or_normal,
equations)
else # u_boundary is "left" of boundary, u_inner is "right" of boundary
flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal,
equations)
noncons = nonconservative_flux_function(u_boundary, u_inner,
orientation_or_normal,
equations)
end

return flux
return flux + 0.5f0 * noncons
end

# Calculate 1D flux for a single point
Expand Down
18 changes: 14 additions & 4 deletions src/equations/shallow_water_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -180,8 +180,10 @@ For details see Section 9.2.5 of the book:
"""
@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector,
x, t,
surface_flux_function,
surface_flux_functions,
equations::ShallowWaterEquations2D)
surface_flux_function, nonconservative_flux_function = surface_flux_functions

# normalize the outward pointing direction
normal = normal_direction / norm(normal_direction)

Expand All @@ -196,8 +198,10 @@ For details see Section 9.2.5 of the book:

# calculate the boundary flux
flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)
noncons = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
equations)

return flux
return flux + 0.5f0 * noncons
end

"""
Expand All @@ -208,8 +212,10 @@ Should be used together with [`TreeMesh`](@ref).
"""
@inline function boundary_condition_slip_wall(u_inner, orientation,
direction, x, t,
surface_flux_function,
surface_flux_functions,
equations::ShallowWaterEquations2D)
# The boundary conditions for the non-conservative term are identically 0 here.
surface_flux_function, nonconservative_flux_function = surface_flux_functions
## get the appropriate normal vector from the orientation
if orientation == 1
u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4])
Expand All @@ -220,11 +226,15 @@ Should be used together with [`TreeMesh`](@ref).
# Calculate boundary flux
if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
noncons = nonconservative_flux_function(u_inner, u_boundary, orientation,
equations)
else # u_boundary is "left" of boundary, u_inner is "right" of boundary
flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
noncons = nonconservative_flux_function(u_boundary, u_inner, orientation,
equations)
end

return flux
return flux + 0.5f0 * noncons
end

# Calculate 1D flux for a single point
Expand Down
17 changes: 4 additions & 13 deletions src/solvers/dgmulti/dg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -519,7 +519,6 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
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.
Expand Down Expand Up @@ -549,18 +548,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)

# Compute pointwise nonconservative numerical flux at the boundary.
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]
dg.surface_integral.surface_flux,
equations)

flux_face_values[i, f] = (cons_flux_at_face_node) * Jf[i, f]
end
end

Expand Down
12 changes: 3 additions & 9 deletions src/solvers/dgsem_p4est/dg_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -349,7 +349,6 @@ end
boundary_index)
@unpack boundaries = cache
@unpack node_coordinates, contravariant_vectors = cache.elements
surface_flux, nonconservative_flux = surface_integral.surface_flux

# Extract solution data from boundary container
u_inner = get_node_vars(boundaries.u, equations, dg, node_index, boundary_index)
Expand All @@ -364,20 +363,15 @@ end

# Call pointwise numerical flux function for the conservative part
# in the normal direction on the boundary
flux_ = boundary_condition(u_inner, normal_direction, x, t, surface_flux, equations)

# Compute pointwise nonconservative numerical flux at the boundary.
noncons_ = boundary_condition(u_inner, normal_direction, x, t, nonconservative_flux,
equations)
flux_ = boundary_condition(u_inner, normal_direction, x, t,
surface_integral.surface_flux, equations)

# Copy flux to element storage in the correct orientation
for v in eachvariable(equations)
# Note the factor 0.5 necessary for the nonconservative fluxes based on
# the interpretation of global SBP operators coupled discontinuously via
# central fluxes/SATs
Comment on lines 371 to 373
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This comment would possibly need added to where the factor of 0.5 scaling occurs.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am also not very content that we have to expose this to the user in order to specify the boundary condition, but I don't see a better way to accomplish this. I definitely think that we should comment this to give some explanation where the factor 0.5 comes from and that it is not specific to any boundary condition.

Copy link
Author

@MarcoArtiano MarcoArtiano Dec 18, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with all of you. I'm also not a huge fan of the factor 0.5. We may alternatively return a tuple from the boundary conditions for non-conservative systems and multiply by 0.5 inside the solver. As an example:

flux_, noncons_ = boundary_condition(u_inner, normal_direction, x, t, surface_integral.surface_flux, equations)

  # Copy flux to element storage in the correct orientation
  for v in eachvariable(equations)
      # Note the factor 0.5 necessary for the nonconservative fluxes based on
      # the interpretation of global SBP operators coupled discontinuously via
      # central fluxes/SATs
      surface_flux_values[v, node_index, direction_index, element_index] = flux_[v] + 0.5f0 * noncons_[v]
  end

What do you think?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I like this idea, since it doesn't require knowledge about the implementation aspect of adding the 0.5 to create the boundary condition. It should work well for the BCs that are currently implemented, but I am not sure what this would look like for the new type of boundary condition that you plan to add. Would you then only set the flux_ part and set the noncons_ part to zero?
I think it would be helpful to include a specific example with such a new BC in the PR to better understand and evaluate what the implementation should look like.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also like this idea with the documentation of how it would work in practice from Patrick's comment above. That is, for something like the standard shallow water equations the jump in the bottom topography is zero at the physical boundary (typically) so one's new boundary condition could compute the conservative flux pieces and then have

   return flux_, SVector(0,0,0,0) # flux, noncons

surface_flux_values[v, node_index, direction_index, element_index] = flux_[v] +
0.5f0 *
noncons_[v]
surface_flux_values[v, node_index, direction_index, element_index] = flux_[v]
end
end

Expand Down
10 changes: 3 additions & 7 deletions src/solvers/dgsem_tree/dg_1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -564,7 +564,6 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
nonconservative_terms::True, equations,
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

@threaded for boundary in first_boundary:last_boundary
Expand All @@ -579,17 +578,14 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
u_inner = u_rr
end
x = get_node_coords(node_coordinates, equations, dg, boundary)

flux = boundary_condition(u_inner, orientations[boundary], direction, x, t,
surface_flux,
surface_integral.surface_flux,
equations)
noncons_flux = boundary_condition(u_inner, orientations[boundary], direction, x,
t, nonconservative_flux,
equations)

# Copy flux to left and right element storage
for v in eachvariable(equations)
surface_flux_values[v, direction, neighbor] = flux[v] +
0.5f0 * noncons_flux[v]
surface_flux_values[v, direction, neighbor] = flux[v]
end
end

Expand Down
9 changes: 2 additions & 7 deletions src/solvers/dgsem_tree/dg_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -762,7 +762,6 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
nonconservative_terms::True, equations,
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

@threaded for boundary in first_boundary:last_boundary
Expand All @@ -779,16 +778,12 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
end
x = get_node_coords(node_coordinates, equations, dg, i, boundary)
flux = boundary_condition(u_inner, orientations[boundary], direction, x, t,
surface_flux,
surface_integral.surface_flux,
equations)
noncons_flux = boundary_condition(u_inner, orientations[boundary],
direction, x, t, nonconservative_flux,
equations)

# Copy flux to left and right element storage
for v in eachvariable(equations)
surface_flux_values[v, i, direction, neighbor] = flux[v] +
0.5f0 * noncons_flux[v]
surface_flux_values[v, i, direction, neighbor] = flux[v]
end
end
end
Expand Down
10 changes: 3 additions & 7 deletions src/solvers/dgsem_unstructured/dg_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -426,7 +426,6 @@ end
surface_integral, dg::DG, cache,
node_index, side_index, element_index,
boundary_index)
surface_flux, nonconservative_flux = surface_integral.surface_flux
@unpack normal_directions = cache.elements
@unpack u, node_coordinates = cache.boundaries

Expand All @@ -442,19 +441,16 @@ end

# Call pointwise numerical flux function for the conservative part
# in the normal direction on the boundary
flux = boundary_condition(u_inner, outward_direction, x, t, surface_flux, equations)
flux = boundary_condition(u_inner, outward_direction, x, t,
surface_integral.surface_flux, equations)

# Compute pointwise nonconservative numerical flux at the boundary.
noncons_flux = boundary_condition(u_inner, outward_direction, x, t,
nonconservative_flux, equations)

for v in eachvariable(equations)
# Note the factor 0.5 necessary for the nonconservative fluxes based on
# the interpretation of global SBP operators coupled discontinuously via
# central fluxes/SATs
surface_flux_values[v, node_index, side_index, element_index] = flux[v] +
0.5f0 *
noncons_flux[v]
surface_flux_values[v, node_index, side_index, element_index] = flux[v]
end
end

Expand Down
4 changes: 2 additions & 2 deletions test/test_type.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2150,7 +2150,7 @@ isdir(outdir) && rm(outdir, recursive = true)
directions = [1, 2]
normal_direction = SVector(one(RealT))

surface_flux_function = flux_lax_friedrichs
surface_flux_function = (flux_lax_friedrichs, flux_lax_friedrichs)
dissipation = DissipationLocalLaxFriedrichs()
numflux = FluxHLL()

Expand Down Expand Up @@ -2239,7 +2239,7 @@ isdir(outdir) && rm(outdir, recursive = true)
directions = [1, 2, 3, 4]
normal_direction = SVector(one(RealT), zero(RealT))

surface_flux_function = flux_lax_friedrichs
surface_flux_function = (flux_lax_friedrichs, flux_lax_friedrichs)
dissipation = DissipationLocalLaxFriedrichs()
numflux = FluxHLL()

Expand Down
Loading