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

Add flexibility to nonconservative BCs #2200

Open
wants to merge 28 commits into
base: main
Choose a base branch
from
+134 −55
Open

Add flexibility to nonconservative BCs #2200

Show file tree
Hide file tree
Changes from 26 commits
Commits
Show all changes
28 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
b60eaf1
moving 0.5 factor internally
MarcoArtiano Jan 7, 2025
d7dce5e
Merge branch 'main' into noncons
MarcoArtiano Jan 8, 2025
840fa9e
Update src/solvers/dgsem_p4est/dg_2d.jl
MarcoArtiano Jan 8, 2025
2baa52a
Update src/equations/shallow_water_2d.jl
MarcoArtiano Jan 8, 2025
cd2186b
Update src/equations/shallow_water_2d.jl
MarcoArtiano Jan 8, 2025
7a326fe
Update src/equations/equations.jl
MarcoArtiano Jan 8, 2025
f9b19c3
Update examples/dgmulti_2d/elixir_mhd_reflective_wall.jl
MarcoArtiano Jan 8, 2025
c35cff9
Update src/equations/shallow_water_2d.jl
MarcoArtiano Jan 8, 2025
26a756a
Update src/equations/shallow_water_2d.jl
MarcoArtiano Jan 8, 2025
307197d
Update src/equations/shallow_water_2d.jl
MarcoArtiano Jan 8, 2025
8a1a3d5
Update src/equations/equations.jl
MarcoArtiano Jan 8, 2025
5068717
Update src/equations/equations.jl
MarcoArtiano Jan 8, 2025
06b881e
Update src/equations/equations.jl
MarcoArtiano Jan 8, 2025
4d112a4
Update src/equations/shallow_water_1d.jl
MarcoArtiano Jan 8, 2025
ecaf306
Update src/equations/shallow_water_1d.jl
MarcoArtiano Jan 8, 2025
1bd0184
Update src/equations/shallow_water_1d.jl
MarcoArtiano Jan 8, 2025
51505aa
minor fixes
MarcoArtiano Jan 8, 2025
d26925e
fix test type stability
MarcoArtiano Jan 8, 2025
9b4e384
fix test type stability
MarcoArtiano Jan 9, 2025
8332bbe
Update src/equations/shallow_water_1d.jl
MarcoArtiano Jan 24, 2025
53d8a14
Update src/equations/shallow_water_2d.jl
MarcoArtiano Jan 24, 2025
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
10 changes: 6 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,10 @@ 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_flux = nonconservative_flux_function(u_inner, u_mirror, normal, equations) *
norm_
return flux, noncons_flux
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)
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
return surface_flux(u_inner, u_inner, orientation_or_normal_direction, equations)
return flux(u_inner, orientation_or_normal_direction, equations)

In a previous PR (https://github.com/trixi-framework/Trixi.jl/pull/2062/files#r1768020588) I changed this from flux to surface_flux. such that the BC works for both conservative and nonconservative systems. Now that we treat them separately, we should be able to only use the flux-function in the conservative case.

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),
Comment on lines +84 to +85
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
return surface_flux_function(u_inner, u_inner, orientation_or_normal_direction,
equations),
return surface_flux_function(u_inner, u_inner,
orientation_or_normal_direction, equations),

For nicer/unified formatting?

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)
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
return surface_flux(u_inner, u_inner, outward_direction, equations)
return flux(u_inner, outward_direction, equations)

see comment above

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),
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 @@ -181,6 +181,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_flux = 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_flux = nonconservative_flux_function(u_boundary, u_inner,
orientation_or_normal,
equations)
end

return flux, noncons_flux
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 @@ -197,6 +226,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_flux = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
equations)
return flux, noncons_flux
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.
patrickersing marked this conversation as resolved.
Show resolved Hide resolved
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_flux = 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_flux = nonconservative_flux_function(u_boundary, u_inner,
orientation_or_normal,
equations)
end

return flux
return flux, noncons_flux
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_flux = nonconservative_flux_function(u_inner, u_boundary, normal_direction,
equations)

return flux
return flux, noncons_flux
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.
MarcoArtiano marked this conversation as resolved.
Show resolved Hide resolved
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_flux = 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_flux = nonconservative_flux_function(u_boundary, u_inner, orientation,
equations)
end

return flux
return flux, noncons_flux
end

# Calculate 1D flux for a single point
Expand Down
22 changes: 8 additions & 14 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 All @@ -546,21 +545,16 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,

# 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.
noncons_flux_at_face_node = boundary_condition(u_face_values[i, f],
face_normal,
face_coordinates,
t,
nonconservative_flux,
equations)
cons_flux_at_face_node, noncons_flux_at_face_node = boundary_condition(u_face_values[i,
f],
face_normal,
face_coordinates,
t,
dg.surface_integral.surface_flux,
equations)

flux_face_values[i, f] = (cons_flux_at_face_node +
0.5 * noncons_flux_at_face_node) * Jf[i, f]
0.5f0 * noncons_flux_at_face_node) * Jf[i, f]
end
end

Expand Down
12 changes: 4 additions & 8 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,17 @@ end

# Call pointwise numerical flux function for the conservative part
# in the normal direction on the boundary
Comment on lines 364 to 365
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
# Call pointwise numerical flux function for the conservative part
# in the normal direction on the boundary
# Call pointwise numerical flux functions for the conservative and nonconservative 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, noncons_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)
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
# 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] +
surface_flux_values[v, node_index, direction_index, element_index] = flux[v] +
0.5f0 *
noncons_[v]
noncons_flux[v]
end
end

Expand Down
12 changes: 5 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,12 +578,11 @@ 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,
equations)
noncons_flux = boundary_condition(u_inner, orientations[boundary], direction, x,
t, nonconservative_flux,
equations)

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

# Copy flux to left and right element storage
for v in eachvariable(equations)
Expand Down
11 changes: 4 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 @@ -778,12 +777,10 @@ function calc_boundary_flux_by_direction!(surface_flux_values::AbstractArray{<:A
u_inner = u_rr
end
x = get_node_coords(node_coordinates, equations, dg, i, boundary)
flux = boundary_condition(u_inner, orientations[boundary], direction, x, t,
surface_flux,
equations)
noncons_flux = boundary_condition(u_inner, orientations[boundary],
direction, x, t, nonconservative_flux,
equations)
flux, noncons_flux = boundary_condition(u_inner, orientations[boundary],
direction, x, t,
surface_integral.surface_flux,
equations)

# Copy flux to left and right element storage
for v in eachvariable(equations)
Expand Down
6 changes: 2 additions & 4 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,11 +441,10 @@ end

# Call pointwise numerical flux function for the conservative part
# in the normal direction on the boundary
Comment on lines 442 to 443
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
# Call pointwise numerical flux function for the conservative part
# in the normal direction on the boundary
# Call pointwise numerical flux functions for the conservative and nonconservative part
# in the normal direction on the boundary

flux = boundary_condition(u_inner, outward_direction, x, t, surface_flux, equations)
flux, noncons_flux = boundary_condition(u_inner, outward_direction, x, t,
surface_integral.surface_flux, equations)

# Compute pointwise nonconservative numerical flux at the boundary.
Comment on lines 446 to 447
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
# Compute pointwise nonconservative numerical flux at the boundary.

This comment can be removed.

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
Expand Down
Loading
Loading