Skip to content

Commit

Permalink
Merge branch 'main' into PERK_special_case
Browse files Browse the repository at this point in the history
  • Loading branch information
warisa-r authored Nov 25, 2024
2 parents f243157 + 62ff764 commit 78df3cc
Show file tree
Hide file tree
Showing 8 changed files with 76 additions and 79 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Trixi"
uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
version = "0.9.7-DEV"
version = "0.9.8-DEV"

[deps]
Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
Expand Down
70 changes: 22 additions & 48 deletions src/callbacks_step/analysis_dg2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -279,30 +279,17 @@ end

function analyze(::Val{:l2_divb}, du, u, t,
mesh::TreeMesh{2},
equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
equations, dg::DGSEM, cache)
integrate_via_indices(u, mesh, equations, dg, cache, cache,
dg.basis.derivative_matrix) do u, i, j, element, equations,
dg, cache, derivative_matrix
divb = zero(eltype(u))
for k in eachnode(dg)
divb += (derivative_matrix[i, k] * u[6, k, j, element] +
derivative_matrix[j, k] * u[7, i, k, element])
end
divb *= cache.elements.inverse_jacobian[element]
divb^2
end |> sqrt
end
B1_kj, _, _ = magnetic_field(u[:, k, j, element], equations)
_, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)

function analyze(::Val{:l2_divb}, du, u, t,
mesh::TreeMesh{2}, equations::IdealGlmMhdMulticomponentEquations2D,
dg::DG, cache)
integrate_via_indices(u, mesh, equations, dg, cache, cache,
dg.basis.derivative_matrix) do u, i, j, element, equations,
dg, cache, derivative_matrix
divb = zero(eltype(u))
for k in eachnode(dg)
divb += (derivative_matrix[i, k] * u[5, k, j, element] +
derivative_matrix[j, k] * u[6, i, k, element])
divb += (derivative_matrix[i, k] * B1_kj +
derivative_matrix[j, k] * B2_ik)
end
divb *= cache.elements.inverse_jacobian[element]
divb^2
Expand All @@ -312,7 +299,7 @@ end
function analyze(::Val{:l2_divb}, du, u, t,
mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
T8codeMesh{2}},
equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
equations, dg::DGSEM, cache)
@unpack contravariant_vectors = cache.elements
integrate_via_indices(u, mesh, equations, dg, cache, cache,
dg.basis.derivative_matrix) do u, i, j, element, equations,
Expand All @@ -323,10 +310,13 @@ function analyze(::Val{:l2_divb}, du, u, t,
Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
# Compute the transformed divergence
for k in eachnode(dg)
B1_kj, B2_kj, _ = magnetic_field(u[:, k, j, element], equations)
B1_ik, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)

divb += (derivative_matrix[i, k] *
(Ja11 * u[6, k, j, element] + Ja12 * u[7, k, j, element]) +
(Ja11 * B1_kj + Ja12 * B2_kj) +
derivative_matrix[j, k] *
(Ja21 * u[6, i, k, element] + Ja22 * u[7, i, k, element]))
(Ja21 * B1_ik + Ja22 * B2_ik))
end
divb *= cache.elements.inverse_jacobian[i, j, element]
divb^2
Expand All @@ -335,7 +325,7 @@ end

function analyze(::Val{:linf_divb}, du, u, t,
mesh::TreeMesh{2},
equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
equations, dg::DGSEM, cache)
@unpack derivative_matrix, weights = dg.basis

# integrate over all elements to get the divergence-free condition errors
Expand All @@ -344,30 +334,11 @@ function analyze(::Val{:linf_divb}, du, u, t,
for j in eachnode(dg), i in eachnode(dg)
divb = zero(eltype(u))
for k in eachnode(dg)
divb += (derivative_matrix[i, k] * u[6, k, j, element] +
derivative_matrix[j, k] * u[7, i, k, element])
end
divb *= cache.elements.inverse_jacobian[element]
linf_divb = max(linf_divb, abs(divb))
end
end

return linf_divb
end

function analyze(::Val{:linf_divb}, du, u, t,
mesh::TreeMesh{2}, equations::IdealGlmMhdMulticomponentEquations2D,
dg::DG, cache)
@unpack derivative_matrix, weights = dg.basis
B1_kj, _, _ = magnetic_field(u[:, k, j, element], equations)
_, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)

# integrate over all elements to get the divergence-free condition errors
linf_divb = zero(eltype(u))
for element in eachelement(dg, cache)
for j in eachnode(dg), i in eachnode(dg)
divb = zero(eltype(u))
for k in eachnode(dg)
divb += (derivative_matrix[i, k] * u[5, k, j, element] +
derivative_matrix[j, k] * u[6, i, k, element])
divb += (derivative_matrix[i, k] * B1_kj +
derivative_matrix[j, k] * B2_ik)
end
divb *= cache.elements.inverse_jacobian[element]
linf_divb = max(linf_divb, abs(divb))
Expand All @@ -380,7 +351,7 @@ end
function analyze(::Val{:linf_divb}, du, u, t,
mesh::Union{StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2},
T8codeMesh{2}},
equations::IdealGlmMhdEquations2D, dg::DGSEM, cache)
equations, dg::DGSEM, cache)
@unpack derivative_matrix, weights = dg.basis
@unpack contravariant_vectors = cache.elements

Expand All @@ -396,10 +367,13 @@ function analyze(::Val{:linf_divb}, du, u, t,
element)
# Compute the transformed divergence
for k in eachnode(dg)
B1_kj, B2_kj, _ = magnetic_field(u[:, k, j, element], equations)
B1_ik, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)

divb += (derivative_matrix[i, k] *
(Ja11 * u[6, k, j, element] + Ja12 * u[7, k, j, element]) +
(Ja11 * B1_kj + Ja12 * B2_kj) +
derivative_matrix[j, k] *
(Ja21 * u[6, i, k, element] + Ja22 * u[7, i, k, element]))
(Ja21 * B1_ik + Ja22 * B2_ik))
end
divb *= cache.elements.inverse_jacobian[i, j, element]
linf_divb = max(linf_divb, abs(divb))
Expand Down
57 changes: 35 additions & 22 deletions src/callbacks_step/analysis_dg3d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -307,16 +307,20 @@ function analyze(::typeof(entropy_timederivative), du, u, t,
end

function analyze(::Val{:l2_divb}, du, u, t,
mesh::TreeMesh{3}, equations::IdealGlmMhdEquations3D,
mesh::TreeMesh{3}, equations,
dg::DGSEM, cache)
integrate_via_indices(u, mesh, equations, dg, cache, cache,
dg.basis.derivative_matrix) do u, i, j, k, element, equations,
dg, cache, derivative_matrix
divb = zero(eltype(u))
for l in eachnode(dg)
divb += (derivative_matrix[i, l] * u[6, l, j, k, element] +
derivative_matrix[j, l] * u[7, i, l, k, element] +
derivative_matrix[k, l] * u[8, i, j, l, element])
B_ljk = magnetic_field(u[:, l, j, k, element], equations)
B_ilk = magnetic_field(u[:, i, l, k, element], equations)
B_ijl = magnetic_field(u[:, i, j, l, element], equations)

divb += (derivative_matrix[i, l] * B_ljk[1] +
derivative_matrix[j, l] * B_ilk[2] +
derivative_matrix[k, l] * B_ijl[3])
end
divb *= cache.elements.inverse_jacobian[element]
divb^2
Expand All @@ -325,7 +329,7 @@ end

function analyze(::Val{:l2_divb}, du, u, t,
mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
equations::IdealGlmMhdEquations3D,
equations,
dg::DGSEM, cache)
@unpack contravariant_vectors = cache.elements
integrate_via_indices(u, mesh, equations, dg, cache, cache,
Expand All @@ -341,23 +345,24 @@ function analyze(::Val{:l2_divb}, du, u, t,
element)
# Compute the transformed divergence
for l in eachnode(dg)
B_ljk = magnetic_field(u[:, l, j, k, element], equations)
B_ilk = magnetic_field(u[:, i, l, k, element], equations)
B_ijl = magnetic_field(u[:, i, j, l, element], equations)

divb += (derivative_matrix[i, l] *
(Ja11 * u[6, l, j, k, element] + Ja12 * u[7, l, j, k, element] +
Ja13 * u[8, l, j, k, element]) +
(Ja11 * B_ljk[1] + Ja12 * B_ljk[2] + Ja13 * B_ljk[3]) +
derivative_matrix[j, l] *
(Ja21 * u[6, i, l, k, element] + Ja22 * u[7, i, l, k, element] +
Ja23 * u[8, i, l, k, element]) +
(Ja21 * B_ilk[1] + Ja22 * B_ilk[2] + Ja23 * B_ilk[3]) +
derivative_matrix[k, l] *
(Ja31 * u[6, i, j, l, element] + Ja32 * u[7, i, j, l, element] +
Ja33 * u[8, i, j, l, element]))
(Ja31 * B_ijl[1] + Ja32 * B_ijl[2] + Ja33 * B_ijl[3]))
end
divb *= cache.elements.inverse_jacobian[i, j, k, element]
divb^2
end |> sqrt
end

function analyze(::Val{:linf_divb}, du, u, t,
mesh::TreeMesh{3}, equations::IdealGlmMhdEquations3D,
mesh::TreeMesh{3}, equations,
dg::DGSEM, cache)
@unpack derivative_matrix, weights = dg.basis

Expand All @@ -367,9 +372,13 @@ function analyze(::Val{:linf_divb}, du, u, t,
for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
divb = zero(eltype(u))
for l in eachnode(dg)
divb += (derivative_matrix[i, l] * u[6, l, j, k, element] +
derivative_matrix[j, l] * u[7, i, l, k, element] +
derivative_matrix[k, l] * u[8, i, j, l, element])
B_ljk = magnetic_field(u[:, l, j, k, element], equations)
B_ilk = magnetic_field(u[:, i, l, k, element], equations)
B_ijl = magnetic_field(u[:, i, j, l, element], equations)

divb += (derivative_matrix[i, l] * B_ljk[1] +
derivative_matrix[j, l] * B_ilk[2] +
derivative_matrix[k, l] * B_ijl[3])
end
divb *= cache.elements.inverse_jacobian[element]
linf_divb = max(linf_divb, abs(divb))
Expand All @@ -386,7 +395,7 @@ end

function analyze(::Val{:linf_divb}, du, u, t,
mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}},
equations::IdealGlmMhdEquations3D,
equations,
dg::DGSEM, cache)
@unpack derivative_matrix, weights = dg.basis
@unpack contravariant_vectors = cache.elements
Expand All @@ -405,12 +414,16 @@ function analyze(::Val{:linf_divb}, du, u, t,
k, element)
# Compute the transformed divergence
for l in eachnode(dg)
divb += (derivative_matrix[i, l] * (Ja11 * u[6, l, j, k, element] +
Ja12 * u[7, l, j, k, element] + Ja13 * u[8, l, j, k, element]) +
derivative_matrix[j, l] * (Ja21 * u[6, i, l, k, element] +
Ja22 * u[7, i, l, k, element] + Ja23 * u[8, i, l, k, element]) +
derivative_matrix[k, l] * (Ja31 * u[6, i, j, l, element] +
Ja32 * u[7, i, j, l, element] + Ja33 * u[8, i, j, l, element]))
B_ljk = magnetic_field(u[:, l, j, k, element], equations)
B_ilk = magnetic_field(u[:, i, l, k, element], equations)
B_ijl = magnetic_field(u[:, i, j, l, element], equations)

divb += (derivative_matrix[i, l] * (Ja11 * B_ljk[1] +
Ja12 * B_ljk[2] + Ja13 * B_ljk[3]) +
derivative_matrix[j, l] * (Ja21 * B_ilk[1] +
Ja22 * B_ilk[2] + Ja23 * B_ilk[3]) +
derivative_matrix[k, l] * (Ja31 * B_ijl[1] +
Ja32 * B_ijl[2] + Ja33 * B_ijl[3]))
end
divb *= cache.elements.inverse_jacobian[i, j, k, element]
linf_divb = max(linf_divb, abs(divb))
Expand Down
3 changes: 3 additions & 0 deletions src/equations/ideal_glm_mhd_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,9 @@ function default_analysis_integrals(::IdealGlmMhdEquations2D)
(entropy_timederivative, Val(:l2_divb), Val(:linf_divb))
end

# Helper function to extract the magnetic field vector from the conservative variables
magnetic_field(u, equations::IdealGlmMhdEquations2D) = SVector(u[6], u[7], u[8])

# Set initial conditions at physical location `x` for time `t`
"""
initial_condition_constant(x, t, equations::IdealGlmMhdEquations2D)
Expand Down
3 changes: 3 additions & 0 deletions src/equations/ideal_glm_mhd_3d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,9 @@ function default_analysis_integrals(::IdealGlmMhdEquations3D)
(entropy_timederivative, Val(:l2_divb), Val(:linf_divb))
end

# Helper function to extract the magnetic field vector from the conservative variables
magnetic_field(u, equations::IdealGlmMhdEquations3D) = SVector(u[6], u[7], u[8])

# Set initial conditions at physical location `x` for time `t`
"""
initial_condition_constant(x, t, equations::IdealGlmMhdEquations3D)
Expand Down
4 changes: 4 additions & 0 deletions src/equations/ideal_glm_mhd_multicomponent_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -101,6 +101,10 @@ function default_analysis_integrals(::IdealGlmMhdMulticomponentEquations2D)
(entropy_timederivative, Val(:l2_divb), Val(:linf_divb))
end

# Helper function to extract the magnetic field vector from the conservative variables
magnetic_field(u, equations::IdealGlmMhdMulticomponentEquations2D) = SVector(u[5], u[6],
u[7])

"""
initial_condition_convergence_test(x, t, equations::IdealGlmMhdMulticomponentEquations2D)
Expand Down
12 changes: 6 additions & 6 deletions src/equations/shallow_water_1d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -201,7 +201,7 @@ end
Non-symmetric two-point volume flux discretizing the nonconservative (source) term
that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref).
Gives entropy conservation and well-balancedness on both the volume and surface when combined with
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:
Expand All @@ -210,7 +210,7 @@ Further details are available in the papers:
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
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)
"""
Expand All @@ -237,7 +237,7 @@ This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces
conservation and well-balancedness.
Further details for the original finite volume formulation are available in
- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
- Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011)
Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography
[DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042)
and for curvilinear 2D case in the paper:
Expand Down Expand Up @@ -310,7 +310,7 @@ is nonzero this should only be used as a surface flux otherwise the scheme will
For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref).
Details are available in Eq. (4.1) in the paper:
- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
- Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011)
Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography
[DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042)
"""
Expand Down Expand Up @@ -582,7 +582,7 @@ end
equations::ShallowWaterEquations1D)
Calculate Roe-averaged velocity `v_roe` and wavespeed `c_roe = sqrt{g * h_roe}`
See for instance equation (62) in
See for instance equation (62) in
- Paul A. Ullrich, Christiane Jablonowski, and Bram van Leer (2010)
High-order finite-volume methods for the shallow-water equations on the sphere
[DOI: 10.1016/j.jcp.2010.04.044](https://doi.org/10.1016/j.jcp.2010.04.044)
Expand Down Expand Up @@ -631,7 +631,7 @@ end
end

# Calculate the error for the "lake-at-rest" test case where H = h+b should
# be a constant value over time.
# be a constant value over time.
@inline function lake_at_rest_error(u, equations::ShallowWaterEquations1D)
h, _, b = u

Expand Down
4 changes: 2 additions & 2 deletions src/equations/shallow_water_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -329,7 +329,7 @@ This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces
conservation and well-balancedness.
Further details for the original finite volume formulation are available in
- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
- Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011)
Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography
[DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042)
and for curvilinear 2D case in the paper:
Expand Down Expand Up @@ -491,7 +491,7 @@ is nonzero this should only be used as a surface flux otherwise the scheme will
For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref).
Details are available in Eq. (4.1) in the paper:
- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
- Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011)
Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography
[DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042)
"""
Expand Down

0 comments on commit 78df3cc

Please sign in to comment.