diff --git a/ext/TrixiConvexECOSExt.jl b/ext/TrixiConvexECOSExt.jl index 8251fe3eed9..9b897436dbb 100644 --- a/ext/TrixiConvexECOSExt.jl +++ b/ext/TrixiConvexECOSExt.jl @@ -60,7 +60,11 @@ function stability_polynomials!(pnoms, consistency_order, num_stage_evals, end # For optimization only the maximum is relevant - return maximum(abs(pnoms)) + if consistency_order - num_stage_evals == 0 + return maximum(abs.(pnoms)) # If there is no variable to optimize, we need to use the broadcast operator. + else + return maximum(abs(pnoms)) + end end #= @@ -151,7 +155,11 @@ function Trixi.bisect_stability_polynomial(consistency_order, num_eig_vals, println("Concluded stability polynomial optimization \n") end - gamma_opt = evaluate(gamma) + if consistency_order - num_stage_evals != 0 + gamma_opt = evaluate(gamma) + else + gamma_opt = nothing # If there is no variable to optimize, return gamma_opt as nothing. + end # Catch case S = 3 (only one opt. variable) if isa(gamma_opt, Number) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 4cd92ce7c5f..e089133fa17 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -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 @@ -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, @@ -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 @@ -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 @@ -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)) @@ -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 @@ -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)) diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl index 85880380c61..3c9c3a69f5e 100644 --- a/src/callbacks_step/analysis_dg3d.jl +++ b/src/callbacks_step/analysis_dg3d.jl @@ -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 @@ -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, @@ -341,15 +345,16 @@ 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 @@ -357,7 +362,7 @@ function analyze(::Val{:l2_divb}, du, u, t, 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 @@ -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)) @@ -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 @@ -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)) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index 9383fdbe961..ba362525228 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -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) diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index 8d7d0461398..d8c859b7a8a 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -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) diff --git a/src/equations/ideal_glm_mhd_multicomponent_2d.jl b/src/equations/ideal_glm_mhd_multicomponent_2d.jl index 99fe391301d..509a2507cd9 100644 --- a/src/equations/ideal_glm_mhd_multicomponent_2d.jl +++ b/src/equations/ideal_glm_mhd_multicomponent_2d.jl @@ -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) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 3c218eee9d9..29642470262 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -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: @@ -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) """ @@ -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: @@ -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) """ @@ -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) @@ -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 diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 7ba36d6ec92..510d4560f9b 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -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: @@ -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) """ diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl index 0ff648556c4..2451680a505 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl @@ -49,15 +49,16 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan, dtmax, dteps, eig_vals; verbose) - monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, - num_stages) - num_monomial_coeffs = length(monomial_coeffs) - @assert num_monomial_coeffs == coeffs_max - A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) - - a_matrix[:, 1] -= A - a_matrix[:, 2] = A + if num_stages != consistency_order + monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, + num_stages) + num_monomial_coeffs = length(monomial_coeffs) + @assert num_monomial_coeffs == coeffs_max + A = compute_a_coeffs(num_stages, stage_scaling_factors, monomial_coeffs) + a_matrix[:, 1] -= A + a_matrix[:, 2] = A + end return a_matrix, c, dt_opt end diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 477748b28c5..02bc3eeba36 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl @@ -33,21 +33,22 @@ function compute_PairedExplicitRK3_butcher_tableau(num_stages, tspan, # Initialize the array of our solution a_unknown = zeros(num_stages - 2) + # Calculate coefficients of the stability polynomial in monomial form + consistency_order = 3 + dtmax = tspan[2] - tspan[1] + dteps = 1.0f-9 + + num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose) + + monomial_coeffs, dt_opt = bisect_stability_polynomial(consistency_order, + num_eig_vals, num_stages, + dtmax, dteps, + eig_vals; verbose) + # Special case of e = 3 - if num_stages == 3 + if num_stages == consistency_order a_unknown = [0.25] # Use classic SSPRK33 (Shu-Osher) Butcher Tableau else - # Calculate coefficients of the stability polynomial in monomial form - consistency_order = 3 - dtmax = tspan[2] - tspan[1] - dteps = 1.0f-9 - - num_eig_vals, eig_vals = filter_eig_vals(eig_vals; verbose) - - monomial_coeffs, dt_opt = bisect_stability_polynomial(consistency_order, - num_eig_vals, num_stages, - dtmax, dteps, - eig_vals; verbose) monomial_coeffs = undo_normalization!(monomial_coeffs, consistency_order, num_stages)