From 335c1978142704e91c5154cfa2c3f7cd51f5bc44 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Rueda-Ram=C3=ADrez?= Date: Tue, 27 Aug 2024 17:36:46 +0200 Subject: [PATCH] Implemented subcell limiting for non-conservative systems with the DGSEM structured solver --- src/equations/ideal_glm_mhd_2d.jl | 137 +++++++++++ .../dg_2d_subcell_limiters.jl | 215 ++++++++++++++++++ 2 files changed, 352 insertions(+) diff --git a/src/equations/ideal_glm_mhd_2d.jl b/src/equations/ideal_glm_mhd_2d.jl index ab2a4b066a..7835eb9b70 100644 --- a/src/equations/ideal_glm_mhd_2d.jl +++ b/src/equations/ideal_glm_mhd_2d.jl @@ -472,6 +472,143 @@ This function is used to compute the subcell fluxes in dg_2d_subcell_limiters.jl return f end +@inline function flux_nonconservative_powell_local_symmetric(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::IdealGlmMhdEquations2D) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr + + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v3_ll = rho_v3_ll / rho_ll + v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll + + # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector + # at the same node location) while `B_dot_n_rr` uses the averaged normal + # direction. The reason for this is that `v_dot_n_ll` depends only on the left + # state and multiplies some gradient while `B_dot_n_rr` is used to compute + # the divergence of B. + B1_avg = (B1_ll + B1_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + B2_avg = (B2_ll + B2_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + psi_avg = (psi_ll + psi_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] + B_dot_n_avg = B1_avg * normal_direction_average[1] + + B2_avg * normal_direction_average[2] + + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + f = SVector(0, + B1_ll * B_dot_n_avg, + B2_ll * B_dot_n_avg, + B3_ll * B_dot_n_avg, + v_dot_B_ll * B_dot_n_avg + v_dot_n_ll * psi_ll * psi_avg, + v1_ll * B_dot_n_avg, + v2_ll * B_dot_n_avg, + v3_ll * B_dot_n_avg, + v_dot_n_ll * psi_avg) + + return f +end + +@inline function flux_nonconservative_powell_local_symmetric(u_ll, u_rr, + normal_direction_ll::AbstractVector, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeLocal, + nonconservative_term::Integer) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr + + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v3_ll = rho_v3_ll / rho_ll + v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll + + # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector + # at the same node location) while `B_dot_n_rr` uses the averaged normal + # direction. The reason for this is that `v_dot_n_ll` depends only on the left + # state and multiplies some gradient while `B_dot_n_rr` is used to compute + # the divergence of B. + if nonconservative_term == 1 + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + f = SVector(0, + B1_ll, + B2_ll, + B3_ll, + v_dot_B_ll, + v1_ll, + v2_ll, + v3_ll, + 0) + else #nonconservative_term == 2 + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2] + f = SVector(0, + 0, + 0, + 0, + v_dot_n_ll * psi_ll, + 0, + 0, + 0, + v_dot_n_ll) + end + + return f +end + +@inline function flux_nonconservative_powell_local_symmetric(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::IdealGlmMhdEquations2D, + nonconservative_type::NonConservativeSymmetric, + nonconservative_term::Integer) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr + + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v3_ll = rho_v3_ll / rho_ll + v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll + + # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector + # at the same node location) while `B_dot_n_rr` uses the averaged normal + # direction. The reason for this is that `v_dot_n_ll` depends only on the left + # state and multiplies some gradient while `B_dot_n_rr` is used to compute + # the divergence of B. + if nonconservative_term == 1 + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + B1_avg = (B1_ll + B1_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + B2_avg = (B2_ll + B2_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + B_dot_n_avg = B1_avg * normal_direction_average[1] + + B2_avg * normal_direction_average[2] + + f = SVector(0, + B_dot_n_avg, + B_dot_n_avg, + B_dot_n_avg, + B_dot_n_avg, + B_dot_n_avg, + B_dot_n_avg, + B_dot_n_avg, + 0) + else #nonconservative_term == 2 + psi_avg = (psi_ll + psi_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + f = SVector(0, + 0, + 0, + 0, + psi_avg, + 0, + 0, + 0, + psi_avg) + end + + return f +end + """ flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations2D) diff --git a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl index e10a529545..fe44e18af2 100644 --- a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl @@ -108,4 +108,219 @@ return nothing end + +# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element +# (**with non-conservative terms**). +# +# See also `flux_differencing_kernel!`. +# +# The calculation of the non-conservative staggered "fluxes" requires non-conservative +# terms that can be written as a product of local and a symmetric contributions. See, e.g., +# +# - Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts +# Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf. +# +@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u, + mesh::StructuredMesh{2}, + nonconservative_terms::True, equations, + volume_flux, dg::DGSEM, element, cache) + (; contravariant_vectors) = cache.elements + (; weights, derivative_split) = dg.basis + (; flux_temp_threaded, flux_nonconservative_temp_threaded) = cache + (; fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded) = cache + + volume_flux_cons, volume_flux_noncons = volume_flux + + flux_temp = flux_temp_threaded[Threads.threadid()] + flux_noncons_temp = flux_nonconservative_temp_threaded[Threads.threadid()] + + fhat_temp = fhat_temp_threaded[Threads.threadid()] + fhat_noncons_temp = fhat_nonconservative_temp_threaded[Threads.threadid()] + phi = phi_threaded[Threads.threadid()] + + # The FV-form fluxes are calculated in a recursive manner, i.e.: + # fhat_(0,1) = w_0 * FVol_0, + # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j, for j=1,...,N-1, + # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i). + + # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated + # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing` + # and saved in in `flux_temp`. + + # Split form volume flux in orientation 1: x direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + # pull the contravariant vectors in each coordinate direction + Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, element) # x direction + + # All diagonal entries of `derivative_split` are zero. Thus, we can skip + # the computation of the diagonal terms. In addition, we use the symmetry + # of the `volume_flux` to save half of the possible two-point flux + # computations. + + # x direction + for ii in (i + 1):nnodes(dg) + u_node_ii = get_node_vars(u, equations, dg, ii, j, element) + # pull the contravariant vectors and compute the average + Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, + element) + Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) + + # compute the contravariant sharp flux in the direction of the averaged contravariant vector + fluxtilde1 = volume_flux(u_node, u_node_ii, Ja1_avg, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], fluxtilde1, + equations, dg, i, j) + multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], fluxtilde1, + equations, dg, ii, j) + for noncons in 1:n_nonconservative_terms(equations) + # We multiply by 0.5 because that is done in other parts of Trixi + fluxtilde1_noncons = volume_flux_noncons(u_node, u_node_ii, Ja1_node, + Ja1_avg, equations, + NonConservativeSymmetric(), + noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5f0 * derivative_split[i, ii], + fluxtilde1_noncons, + equations, dg, noncons, i, j) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5f0 * derivative_split[ii, i], + fluxtilde1_noncons, + equations, dg, noncons, ii, j) + end + end + end + + # FV-form flux `fhat` in x direction + fhat1_L[:, 1, :] .= zero(eltype(fhat1_L)) + fhat1_L[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_L)) + fhat1_R[:, 1, :] .= zero(eltype(fhat1_R)) + fhat1_R[:, nnodes(dg) + 1, :] .= zero(eltype(fhat1_R)) + + fhat_temp[:, 1, :] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, 1, :] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, element) # x direction + for noncons in 1:n_nonconservative_terms(equations) + set_node_vars!(phi, + volume_flux_noncons(u_local, Ja1_node, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j) + end + end + + for j in eachnode(dg), i in 1:(nnodes(dg) - 1) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j] + weights[i] * flux_temp[v, i, j] + fhat_temp[v, i + 1, j] = value + fhat1_L[v, i + 1, j] = value + fhat1_R[v, i + 1, j] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(equations), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j] + + weights[i] * flux_noncons_temp[v, noncons, i, j] + fhat_noncons_temp[v, noncons, i + 1, j] = value + + fhat1_L[v, i + 1, j] = fhat1_L[v, i + 1, j] + phi[v, noncons, i, j] * value + fhat1_R[v, i + 1, j] = fhat1_R[v, i + 1, j] + + phi[v, noncons, i + 1, j] * value + end + end + + # Split form volume flux in orientation 2: y direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + # pull the contravariant vectors in each coordinate direction + Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element) + + # y direction + for jj in (j + 1):nnodes(dg) + u_node_jj = get_node_vars(u, equations, dg, i, jj, element) + # pull the contravariant vectors and compute the average + Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, + element) + Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) + # compute the contravariant sharp flux in the direction of the averaged contravariant vector + fluxtilde2 = volume_flux(u_node, u_node_jj, Ja2_avg, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], fluxtilde2, + equations, dg, i, j) + multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], fluxtilde2, + equations, dg, i, jj) + for noncons in 1:n_nonconservative_terms(equations) + # We multiply by 0.5 because that is done in other parts of Trixi + fluxtilde2_noncons = volume_flux_noncons(u_node, u_node_jj, Ja2_node, + Ja2_avg, equations, + NonConservativeSymmetric(), + noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5f0 * derivative_split[j, jj], + fluxtilde2_noncons, + equations, dg, noncons, i, j) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5f0 * derivative_split[jj, j], + fluxtilde2_noncons, + equations, dg, noncons, i, jj) + end + end + end + + # FV-form flux `fhat` in y direction + fhat2_L[:, :, 1] .= zero(eltype(fhat2_L)) + fhat2_L[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2_L)) + fhat2_R[:, :, 1] .= zero(eltype(fhat2_R)) + fhat2_R[:, :, nnodes(dg) + 1] .= zero(eltype(fhat2_R)) + + fhat_temp[:, 1, :] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, 1, :] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element) + for noncons in 1:n_nonconservative_terms(equations) + set_node_vars!(phi, + volume_flux_noncons(u_local, Ja2_node, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j) + end + end + + for j in 1:(nnodes(dg) - 1), i in eachnode(dg) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j] + weights[j] * flux_temp[v, i, j] + fhat_temp[v, i, j + 1] = value + fhat1_L[v, i, j + 1] = value + fhat1_R[v, i, j + 1] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(equations), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j] + + weights[j] * flux_noncons_temp[v, noncons, i, j] + fhat_noncons_temp[v, noncons, i, j + 1] = value + + fhat1_L[v, i, j + 1] = fhat1_L[v, i, j + 1] + phi[v, noncons, i, j] * value + fhat1_R[v, i, j + 1] = fhat1_R[v, i, j + 1] + + phi[v, noncons, i, j + 1] * value + end + end + + return nothing +end end # @muladd