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

Update NC-fluxes for the SWE #2038

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
4 changes: 2 additions & 2 deletions examples/tree_2d_dgsem/elixir_shallowwater_wall.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,8 @@ boundary_condition = boundary_condition_slip_wall
###############################################################################
# Get the DG approximation space

volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal)
volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
surface_flux = (flux_lax_friedrichs, flux_nonconservative_wintermeyer_etal)
solver = DGSEM(polydeg = 3, surface_flux = surface_flux,
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

Expand Down
1 change: 0 additions & 1 deletion src/Trixi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -179,7 +179,6 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
flux_kennedy_gruber, flux_shima_etal, flux_ec,
flux_fjordholm_etal, flux_nonconservative_fjordholm_etal,
flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal,
flux_nonconservative_ersing_etal,
flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal,
hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs,
Expand Down
69 changes: 16 additions & 53 deletions src/equations/shallow_water_1d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -201,20 +201,27 @@ end
Non-symmetric two-point volume flux discretizing the nonconservative (source) term
that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref).

Further details are available in the paper:
Gives entropy conservation and well-balancedness on both the volume and surface when combined with
[`flux_wintermeyer_etal`](@ref).

Further details are available in the papers:
- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
An entropy stable nodal discontinuous Galerkin method for the two dimensional
shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
[DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
- Patrick Ersing, Andrew R. Winters (2023)
An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on
curvilinear meshes
[DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
"""
@inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
equations::ShallowWaterEquations1D)
# Pull the necessary left and right state information
h_ll = waterheight(u_ll, equations)
b_rr = u_rr[3]
b_jump = u_rr[3] - u_ll[3]

# Bottom gradient nonconservative term: (0, g h b_x, 0)
f = SVector(0, equations.gravity * h_ll * b_rr, 0)
f = SVector(0, equations.gravity * h_ll * b_jump, 0)

return f
end
Expand All @@ -226,11 +233,8 @@ end
Non-symmetric two-point surface flux discretizing the nonconservative (source) term of
that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref).

This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref)
that account for possible discontinuities in the bottom topography function.
Thus, this flux should be used in general at interfaces. For flux differencing volume terms,
[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly
cheaper.
This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy
conservation and well-balancedness.

Further details for the original finite volume formulation are available in
- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
Expand All @@ -256,7 +260,6 @@ and for curvilinear 2D case in the paper:
# cross-averaging across a discontinuous bottom topography
# (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry
f = SVector(0,
equations.gravity * h_ll * b_ll +
equations.gravity * h_average * b_jump,
0)

Expand Down Expand Up @@ -285,60 +288,19 @@ Further details on the hydrostatic reconstruction and its motivation can be foun
orientation::Integer,
equations::ShallowWaterEquations1D)
# Pull the water height and bottom topography on the left
h_ll, _, b_ll = u_ll
h_ll, _, _ = u_ll

# Create the hydrostatic reconstruction for the left solution state
u_ll_star, _ = hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, equations)

# Copy the reconstructed water height for easier to read code
h_ll_star = u_ll_star[1]

# Includes two parts:
# (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
# cross-averaging across a discontinuous bottom topography
# (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry
return SVector(0,
equations.gravity * h_ll * b_ll +
equations.gravity * (h_ll^2 - h_ll_star^2),
0)
end

"""
flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
equations::ShallowWaterEquations1D)

!!! warning "Experimental code"
This numerical flux is experimental and may change in any future release.

Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term
that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref).

This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy
conservation and well-balancedness in both the volume and surface when combined with
[`flux_wintermeyer_etal`](@ref).

For further details see:
- Patrick Ersing, Andrew R. Winters (2023)
An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on
curvilinear meshes
[DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
"""
@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
equations::ShallowWaterEquations1D)
# Pull the necessary left and right state information
h_ll = waterheight(u_ll, equations)
b_rr = u_rr[3]
b_ll = u_ll[3]

# Calculate jump
b_jump = b_rr - b_ll

# Bottom gradient nonconservative term: (0, g h b_x, 0)
f = SVector(0, equations.gravity * h_ll * b_jump, 0)

return f
end

"""
flux_fjordholm_etal(u_ll, u_rr, orientation,
equations::ShallowWaterEquations1D)
Expand Down Expand Up @@ -378,7 +340,8 @@ end

Total energy conservative (mathematical entropy for shallow water equations) split form.
When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref).
For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can
be used to ensure well-balancedness and entropy conservation.

Further details are available in Theorem 1 of the paper:
- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
Expand Down Expand Up @@ -571,7 +534,7 @@ end
# Note, only the first two are the entropy variables, the third entry still
# just carries the bottom topography values for convenience
@inline function cons2entropy(u, equations::ShallowWaterEquations1D)
h, h_v, b = u
h, _, b = u

v = velocity(u, equations)

Expand Down
151 changes: 25 additions & 126 deletions src/equations/shallow_water_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -274,28 +274,30 @@ end
Non-symmetric two-point volume flux discretizing the nonconservative (source) term
that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).

On curvilinear meshes, this nonconservative flux depends on both the
contravariant vector (normal direction) at the current node and the averaged
one. This is different from numerical fluxes used to discretize conservative
terms.
For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can
be used to ensure well-balancedness and entropy conservation.

Further details are available in the paper:
Further details are available in the papers:
- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
An entropy stable nodal discontinuous Galerkin method for the two dimensional
shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
[DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
- Patrick Ersing, Andrew R. Winters (2023)
An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on
curvilinear meshes
[DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
"""
@inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
equations::ShallowWaterEquations2D)
# Pull the necessary left and right state information
h_ll = waterheight(u_ll, equations)
b_rr = u_rr[4]
b_jump = u_rr[4] - u_ll[4]

# Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
if orientation == 1
f = SVector(0, equations.gravity * h_ll * b_rr, 0, 0)
f = SVector(0, equations.gravity * h_ll * b_jump, 0, 0)
else # orientation == 2
f = SVector(0, 0, equations.gravity * h_ll * b_rr, 0)
f = SVector(0, 0, equations.gravity * h_ll * b_jump, 0)
end
return f
end
Expand All @@ -306,12 +308,12 @@ end
equations::ShallowWaterEquations2D)
# Pull the necessary left and right state information
h_ll = waterheight(u_ll, equations)
b_rr = u_rr[4]
# Note this routine only uses the `normal_direction_average` and the average of the
# bottom topography to get a quadratic split form DG gradient on curved elements
b_jump = u_rr[4] - u_ll[4]

# Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
return SVector(0,
normal_direction_average[1] * equations.gravity * h_ll * b_rr,
normal_direction_average[2] * equations.gravity * h_ll * b_rr,
normal_direction_average[1] * equations.gravity * h_ll * b_jump,
normal_direction_average[2] * equations.gravity * h_ll * b_jump,
0)
end

Expand All @@ -326,16 +328,8 @@ end
Non-symmetric two-point surface flux discretizing the nonconservative (source) term of
that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).

On curvilinear meshes, this nonconservative flux depends on both the
contravariant vector (normal direction) at the current node and the averaged
one. This is different from numerical fluxes used to discretize conservative
terms.

This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref)
that account for possible discontinuities in the bottom topography function.
Thus, this flux should be used in general at interfaces. For flux differencing volume terms,
[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly
cheaper.
This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy
conservation and well-balancedness.

Further details for the original finite volume formulation are available in
- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
Expand All @@ -356,18 +350,13 @@ and for curvilinear 2D case in the paper:
h_average = 0.5f0 * (h_ll + h_rr)
b_jump = b_rr - b_ll

# Includes two parts:
# (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
# cross-averaging across a discontinuous bottom topography
# (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry
# Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
if orientation == 1
f = SVector(0,
equations.gravity * h_ll * b_ll +
equations.gravity * h_average * b_jump,
0, 0)
else # orientation == 2
f = SVector(0, 0,
equations.gravity * h_ll * b_ll +
equations.gravity * h_average * b_jump,
0)
end
Expand All @@ -383,20 +372,12 @@ end
h_ll, _, _, b_ll = u_ll
h_rr, _, _, b_rr = u_rr

# Comes in two parts:
# (i) Diagonal (consistent) term from the volume flux that uses `normal_direction_average`
# but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography

f2 = normal_direction_average[1] * equations.gravity * h_ll * b_ll
f3 = normal_direction_average[2] * equations.gravity * h_ll * b_ll

# (ii) True surface part that uses `normal_direction_ll`, `h_average` and `b_jump`
# to handle discontinuous bathymetry
h_average = 0.5f0 * (h_ll + h_rr)
b_jump = b_rr - b_ll

f2 += normal_direction_ll[1] * equations.gravity * h_average * b_jump
f3 += normal_direction_ll[2] * equations.gravity * h_average * b_jump
# Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
f2 = normal_direction_average[1] * equations.gravity * h_average * b_jump
f3 = normal_direction_average[2] * equations.gravity * h_average * b_jump

# First and last equations do not have a nonconservative flux
f1 = f4 = 0
Expand Down Expand Up @@ -472,18 +453,12 @@ Further details for the hydrostatic reconstruction and its motivation can be fou
# Copy the reconstructed water height for easier to read code
h_ll_star = u_ll_star[1]

# Includes two parts:
# (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
# cross-averaging across a discontinuous bottom topography
# (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry
if orientation == 1
f = SVector(0,
equations.gravity * h_ll * b_ll +
equations.gravity * (h_ll^2 - h_ll_star^2),
0, 0)
else # orientation == 2
f = SVector(0, 0,
equations.gravity * h_ll * b_ll +
equations.gravity * (h_ll^2 - h_ll_star^2),
0)
end
Expand All @@ -504,92 +479,15 @@ end
# Copy the reconstructed water height for easier to read code
h_ll_star = u_ll_star[1]

# Comes in two parts:
# (i) Diagonal (consistent) term from the volume flux that uses `normal_direction_average`
# but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography

f2 = normal_direction_average[1] * equations.gravity * h_ll * b_ll
f3 = normal_direction_average[2] * equations.gravity * h_ll * b_ll

# (ii) True surface part that uses `normal_direction_ll`, `h_ll` and `h_ll_star`
# to handle discontinuous bathymetry

f2 += normal_direction_ll[1] * equations.gravity * (h_ll^2 - h_ll_star^2)
f3 += normal_direction_ll[2] * equations.gravity * (h_ll^2 - h_ll_star^2)
f2 = normal_direction_average[1] * equations.gravity * (h_ll^2 - h_ll_star^2)
f3 = normal_direction_average[2] * equations.gravity * (h_ll^2 - h_ll_star^2)

# First and last equations do not have a nonconservative flux
f1 = f4 = 0

return SVector(f1, f2, f3, f4)
end

"""
flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
equations::ShallowWaterEquations2D)
flux_nonconservative_ersing_etal(u_ll, u_rr,
normal_direction_ll::AbstractVector,
normal_direction_average::AbstractVector,
equations::ShallowWaterEquations2D)

!!! warning "Experimental code"
This numerical flux is experimental and may change in any future release.

Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term
that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).

On curvilinear meshes, this nonconservative flux depends on both the
contravariant vector (normal direction) at the current node and the averaged
one. This is different from numerical fluxes used to discretize conservative
terms.

This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy
conservation and well-balancedness in both the volume and surface when combined with
[`flux_wintermeyer_etal`](@ref).

For further details see:
- Patrick Ersing, Andrew R. Winters (2023)
An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on
curvilinear meshes
[DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
"""
@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
equations::ShallowWaterEquations2D)
# Pull the necessary left and right state information
h_ll = waterheight(u_ll, equations)
b_rr = u_rr[4]
b_ll = u_ll[4]

# Calculate jump
b_jump = b_rr - b_ll

# Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
if orientation == 1
f = SVector(0, equations.gravity * h_ll * b_jump, 0, 0)
else # orientation == 2
f = SVector(0, 0, equations.gravity * h_ll * b_jump, 0)
end
return f
end

@inline function flux_nonconservative_ersing_etal(u_ll, u_rr,
normal_direction_ll::AbstractVector,
normal_direction_average::AbstractVector,
equations::ShallowWaterEquations2D)
# Pull the necessary left and right state information
h_ll = waterheight(u_ll, equations)
b_rr = u_rr[4]
b_ll = u_ll[4]

# Calculate jump
b_jump = b_rr - b_ll
# Note this routine only uses the `normal_direction_average` and the average of the
# bottom topography to get a quadratic split form DG gradient on curved elements
return SVector(0,
normal_direction_average[1] * equations.gravity * h_ll * b_jump,
normal_direction_average[2] * equations.gravity * h_ll * b_jump,
0)
end

"""
flux_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction,
equations::ShallowWaterEquations2D)
Expand Down Expand Up @@ -664,7 +562,8 @@ end

Total energy conservative (mathematical entropy for shallow water equations) split form.
When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref).
For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can
be used to ensure well-balancedness and entropy conservation.

Further details are available in Theorem 1 of the paper:
- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
Expand Down
Loading
Loading