EC mortars: Tree 2D

examples/tree_2d_dgsem/elixir_euler_vortex_mortar_ec.jl

@@ -0,0 +1,106 @@
+
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+"""
+    initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)
+
+The classical isentropic vortex test case of
+- Chi-Wang Shu (1997)
+  Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory
+  Schemes for Hyperbolic Conservation Laws
+  [NASA/CR-97-206253](https://ntrs.nasa.gov/citations/19980007543)
+"""
+function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)
+    # needs appropriate mesh size, e.g. [-10,-10]x[10,10]
+    # for error convergence: make sure that the end time is such that the vortex is back at the initial state!!
+    # for the current velocity and domain size: t_end should be a multiple of 20s
+    # initial center of the vortex
+    inicenter = SVector(0.0, 0.0)
+    # size and strength of the vortex
+    iniamplitude = 5.0
+    # base flow
+    rho = 1.0
+    v1 = 1.0
+    v2 = 1.0
+    vel = SVector(v1, v2)
+    p = 25.0
+    rt = p / rho                  # ideal gas equation
+    t_loc = 0
+    cent = inicenter + vel * t_loc      # advection of center
+    # ATTENTION: handle periodic BC, but only for v1 = v2 = 1.0 (!!!!)
+
+    cent = x - cent # distance to center point
+
+    # cent = cross(iniaxis, cent) # distance to axis, tangent vector, length r
+    # cross product with iniaxis = [0, 0, 1]
+    cent = SVector(-cent[2], cent[1])
+    r2 = cent[1]^2 + cent[2]^2
+    du = iniamplitude / (2 * π) * exp(0.5 * (1 - r2)) # vel. perturbation
+    dtemp = -(equations.gamma - 1) / (2 * equations.gamma * rt) * du^2 # isentropic
+    rho = rho * (1 + dtemp)^(1 / (equations.gamma - 1))
+    vel = vel + du * cent
+    v1, v2 = vel
+    p = p * (1 + dtemp)^(equations.gamma / (equations.gamma - 1))
+    prim = SVector(rho, v1, v2, p)
+    return prim2cons(prim, equations)
+end
+initial_condition = initial_condition_isentropic_vortex
+
+polydeg = 3
+
+#surface_flux = flux_ranocha # For truly entropy-conservative spatial discretization
+surface_flux = flux_lax_friedrichs # For convergence test
+
+basis = LobattoLegendreBasis(Float64, polydeg)
+solver = DGSEM(basis, surface_flux,
+               VolumeIntegralFluxDifferencing(flux_ranocha),
+               MortarEC(basis))
+
+coordinates_min = (-10.0, -10.0)
+coordinates_max = (10.0, 10.0)
+
+# Add refinement patch 
+refinement_patch = ((type = "box", coordinates_min = (-5.0, -5.0),
+                     coordinates_max = (5.0, 5.0)),)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                refinement_patches = refinement_patch,
+                n_cells_max = 100_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 20.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     save_analysis = true,
+                                     analysis_filename = "analysis.dat",
+                                     extra_analysis_integrals = (entropy,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.4)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);

src/Trixi.jl

@@ -250,7 +250,7 @@ export DG,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        SurfaceIntegralUpwind,
-       MortarL2
+       MortarL2, MortarEC
 
 export VolumeIntegralSubcellLimiting, BoundsCheckCallback,
        SubcellLimiterIDP, SubcellLimiterIDPCorrection

src/solvers/dgsem/basis_lobatto_legendre.jl

@@ -155,6 +155,24 @@ struct LobattoLegendreMortarL2{RealT <: Real, NNODES,
     reverse_lower::ReverseMatrix
 end
 
+"""
+    MortarL2(basis::LobattoLegendreBasis)
+
+L2-projection based mortar tailored to the [`LobattoLegendreBasis`](@ref).
+Mortars are required to handle non-conforming interfaces. 
+Currently, Trixi.jl only supports 2h-nonconforming interfaces.
+This means that all nonconforming interfaces (lines in 2D and faces in 3D) are of a 2-1 (or 1-2) coupling, i.e.,
+side A has exactly twice the number of nodes (per dimension) of side B (or vice-versa).
+
+# References
+- David A. Kopriva, (2009). 
+  Implementing spectral methods for partial differential equations:
+  Algorithms for scientists and engineers. 
+  [DOI:10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
+- David A. Kopriva, (1996). 
+  A Conservative Staggered-Grid Chebyshev Multidomain Method for Compressible Flows. II. A Semi-Structured Method.
+  [DOI:10.1006/jcph.1996.0225](https://doi.org/10.1006/jcph.1996.0225)
+"""
 function MortarL2(basis::LobattoLegendreBasis)
     RealT = real(basis)
     nnodes_ = nnodes(basis)
@@ -196,47 +214,82 @@ function Base.show(io::IO, ::MIME"text/plain", mortar::LobattoLegendreMortarL2)
           polydeg(mortar))
 end
 
-@inline Base.real(mortar::LobattoLegendreMortarL2{RealT}) where {RealT} = RealT
+abstract type AbstractMortarEC{RealT} <: AbstractMortar{RealT} end
+
+struct LobattoLegendreMortarEC{RealT <: Real, NNODES,
+                               ForwardMatrix <: AbstractMatrix{RealT},
+                               ReverseMatrix <: AbstractMatrix{RealT}} <:
+       AbstractMortarEC{RealT}
+    forward_upper::ForwardMatrix
+    forward_lower::ForwardMatrix
+    reverse_upper::ReverseMatrix
+    reverse_lower::ReverseMatrix
+end
+
+"""
+    MortarEC(basis::LobattoLegendreBasis)
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in any future releases.
+    Currently only implemented for 2D [`TreeMesh`](@ref).
+
+Flux-corrected, entropy-conservative L2-projection based mortar tailored to the [`LobattoLegendreBasis`](@ref).
+Mortars are required to handle non-conforming interfaces. 
+Currently, Trixi.jl only supports 2h-nonconforming interfaces.
+This means that all nonconforming interfaces (lines in 2D and faces in 3D) are of a 2-1 (or 1-2) coupling, i.e.,
+side A has exactly twice the number of nodes (per dimension) of side B (or vice-versa).
+
+The entropy-conservative mortar is based on the [`MortarL2`](@ref) but adds a flux correction term
+to ensure that the mortar projection/interpolation operations are entropy-conservative.
+
+# References
+- Jesse Chan, Mario J. Bencomo, and David C. Del Rey Fernández (2021).
+  Mortar-based Entropy-Stable Discontinuous Galerkin Methods on Non-conforming Quadrilateral and Hexahedral Meshes.
+  [DOI:10.1007/s10915-021-01652-3](https://doi.org/10.1007/s10915-021-01652-3)
+"""
+function MortarEC(basis::LobattoLegendreBasis)
+    RealT = real(basis)
+    nnodes_ = nnodes(basis)
+
+    forward_upper = calc_forward_upper(nnodes_, RealT)
+    forward_lower = calc_forward_lower(nnodes_, RealT)
+    reverse_upper = calc_reverse_upper(nnodes_, Val(:gauss), RealT)
+    reverse_lower = calc_reverse_lower(nnodes_, Val(:gauss), RealT)
+
+    LobattoLegendreMortarEC{RealT, nnodes_, typeof(forward_upper),
+                            typeof(reverse_upper)}(forward_upper, forward_lower,
+                                                   reverse_upper, reverse_lower)
+end
+
+function Base.show(io::IO, mortar::LobattoLegendreMortarEC)
+    @nospecialize mortar # reduce precompilation time
+
+    print(io, "LobattoLegendreMortarEC{", real(mortar), "}(polydeg=", polydeg(mortar),
+          ")")
+end
+function Base.show(io::IO, ::MIME"text/plain", mortar::LobattoLegendreMortarEC)
+    @nospecialize mortar # reduce precompilation time
+
+    print(io, "LobattoLegendreMortarEC{", real(mortar), "} with polynomials of degree ",
+          polydeg(mortar))
+end
+
+# Shared type alias for both mortar types to enable common functions in a compact way
+const LobattoLegendreMortar{RealT, NNODES} = Union{LobattoLegendreMortarL2{RealT,
+                                                                           NNODES},
+                                                   LobattoLegendreMortarEC{RealT,
+                                                                           NNODES}}
+
+@inline Base.real(mortar::LobattoLegendreMortar{RealT}) where {RealT} = RealT
 
-@inline function nnodes(mortar::LobattoLegendreMortarL2{RealT, NNODES}) where {RealT,
-                                                                               NNODES}
+@inline function nnodes(mortar::LobattoLegendreMortar{RealT, NNODES}) where {
+                                                                             RealT,
+                                                                             NNODES
+                                                                             }
     NNODES
 end
 
-@inline polydeg(mortar::LobattoLegendreMortarL2) = nnodes(mortar) - 1
-
-# TODO: We can create EC mortars along the lines of the following implementation.
-# abstract type AbstractMortarEC{RealT} <: AbstractMortar{RealT} end
-
-# struct LobattoLegendreMortarEC{RealT<:Real, NNODES, MortarMatrix<:AbstractMatrix{RealT}, SurfaceFlux} <: AbstractMortarEC{RealT}
-#   forward_upper::MortarMatrix
-#   forward_lower::MortarMatrix
-#   reverse_upper::MortarMatrix
-#   reverse_lower::MortarMatrix
-#   surface_flux::SurfaceFlux
-# end
-
-# function MortarEC(basis::LobattoLegendreBasis{RealT}, surface_flux)
-#   forward_upper   = calc_forward_upper(n_nodes, RealT)
-#   forward_lower   = calc_forward_lower(n_nodes, RealT)
-#   l2reverse_upper = calc_reverse_upper(n_nodes, Val(:gauss_lobatto), RealT)
-#   l2reverse_lower = calc_reverse_lower(n_nodes, Val(:gauss_lobatto), RealT)
-
-#   # type conversions to make use of StaticArrays etc.
-#   nnodes_ = nnodes(basis)
-#   forward_upper   = SMatrix{nnodes_, nnodes_}(forward_upper)
-#   forward_lower   = SMatrix{nnodes_, nnodes_}(forward_lower)
-#   l2reverse_upper = SMatrix{nnodes_, nnodes_}(l2reverse_upper)
-#   l2reverse_lower = SMatrix{nnodes_, nnodes_}(l2reverse_lower)
-
-#   LobattoLegendreMortarEC{RealT, nnodes_, typeof(forward_upper), typeof(surface_flux)}(
-#     forward_upper, forward_lower,
-#     l2reverse_upper, l2reverse_lower,
-#     surface_flux
-#   )
-# end
-
-# @inline nnodes(mortar::LobattoLegendreMortarEC{RealT, NNODES}) = NNODES
+@inline polydeg(mortar::LobattoLegendreMortar) = nnodes(mortar) - 1
 
 struct LobattoLegendreAnalyzer{RealT <: Real, NNODES,
                                VectorT <: AbstractVector{RealT},

src/solvers/dgsem_p4est/dg_2d.jl

@@ -444,7 +444,7 @@ function prolong2mortars!(cache, u,
     return nothing
 end
 
-function calc_mortar_flux!(surface_flux_values,
+function calc_mortar_flux!(surface_flux_values, u,
                            mesh::Union{P4estMesh{2}, T8codeMesh{2}},
                            nonconservative_terms, equations,
                            mortar_l2::LobattoLegendreMortarL2,

src/solvers/dgsem_p4est/dg_2d_parabolic.jl

@@ -260,7 +260,7 @@ function calc_gradient!(gradients, u_transformed, t,
     # along with a specialization on `calc_mortar_flux!(fstar, ...)` and `mortar_fluxes_to_elements!` for
     # AbstractEquationsParabolic.
     @trixi_timeit timer() "mortar flux" begin
-        calc_mortar_flux!(cache_parabolic.elements.surface_flux_values,
+        calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, u_transformed,
                           mesh, False(), # False() = no nonconservative terms
                           equations_parabolic,
                           dg.mortar, dg.surface_integral, dg, cache)

src/solvers/dgsem_p4est/dg_3d.jl

@@ -519,7 +519,7 @@ function prolong2mortars!(cache, u,
     return nothing
 end
 
-function calc_mortar_flux!(surface_flux_values,
+function calc_mortar_flux!(surface_flux_values, u,
                            mesh::Union{P4estMesh{3}, T8codeMesh{3}},
                            nonconservative_terms, equations,
                            mortar_l2::LobattoLegendreMortarL2,

src/solvers/dgsem_p4est/dg_3d_parabolic.jl

@@ -155,7 +155,7 @@ function calc_gradient!(gradients, u_transformed, t,
     # along with a specialization on `calc_mortar_flux!` and `mortar_fluxes_to_elements!` for
     # AbstractEquationsParabolic.
     @trixi_timeit timer() "mortar flux" begin
-        calc_mortar_flux!(cache_parabolic.elements.surface_flux_values,
+        calc_mortar_flux!(cache_parabolic.elements.surface_flux_values, u_transformed,
                           mesh, False(), # False() = no nonconservative terms
                           equations_parabolic,
                           dg.mortar, dg.surface_integral, dg, cache)

src/solvers/dgsem_p4est/dg_3d_parallel.jl

@@ -69,8 +69,8 @@ function rhs!(du, u, t,
 
     # Calculate mortar fluxes
     @trixi_timeit timer() "mortar flux" begin
-        calc_mortar_flux!(cache.elements.surface_flux_values, mesh,
-                          have_nonconservative_terms(equations), equations,
+        calc_mortar_flux!(cache.elements.surface_flux_values, u,
+                          mesh, have_nonconservative_terms(equations), equations,
                           dg.mortar, dg.surface_integral, dg, cache)
     end
 

src/solvers/dgsem_tree/containers_2d.jl

@@ -600,9 +600,10 @@ end
 # Create mortar container and initialize mortar data in `elements`.
 function init_mortars(cell_ids, mesh::TreeMesh2D,
                       elements::ElementContainer2D,
-                      ::LobattoLegendreMortarL2)
+                      ::LobattoLegendreMortar)
     # Initialize containers
     n_mortars = count_required_mortars(mesh, cell_ids)
+    # We can reuse the `L2MortarContainer2D` also for the `LobattoLegendreMortarEC` mortar
     mortars = L2MortarContainer2D{eltype(elements)}(n_mortars, nvariables(elements),
                                                     nnodes(elements))