tetrahedral metrics use p+1 interpolant

src/Analysis/conservation.jl

@@ -4,34 +4,35 @@ abstract type AbstractConservationAnalysisResults <: AbstractAnalysisResults end
 """
 Evaluate change in ∫udx  
 """
-struct PrimaryConservationAnalysis <: ConservationAnalysis
-    WJ::Vector{Diagonal}
+struct PrimaryConservationAnalysis{V_type} <: ConservationAnalysis
+    WJ::Vector{Diagonal{Float64, Vector{Float64}}}
     N_c::Int
     N_e::Int
-    V::LinearMap
+    V::V_type
     results_path::String
     dict_name::String
 end
 
 """
 Evaluate change in  ∫½u²dx  
 """
-struct EnergyConservationAnalysis <: ConservationAnalysis
-    mass_solver::AbstractMassMatrixSolver
+struct EnergyConservationAnalysis{V_type, MassSolver} <: ConservationAnalysis
+    mass_solver::MassSolver
     N_c::Int
     N_e::Int
-    V::LinearMap
+    V::V_type
     results_path::String
     dict_name::String
 end
 
-struct EntropyConservationAnalysis <: ConservationAnalysis
-    mass_solver::AbstractMassMatrixSolver
-    WJ::Vector{Diagonal}
-    conservation_law::AbstractConservationLaw
+struct EntropyConservationAnalysis{V_type, MassSolver, 
+    ConservationLaw} <: ConservationAnalysis
+    mass_solver::MassSolver
+    WJ::Vector{Diagonal{Float64, Vector{Float64}}}
+    conservation_law::ConservationLaw
     N_c::Int
     N_e::Int
-    V::LinearMap
+    V::V_type
     results_path::String
     dict_name::String
 end
@@ -285,20 +286,24 @@ end
 end
 
 @recipe function plot(results::ConservationAnalysisResultsWithDerivative,
-    e::Int=1)
+    e::Int=1; net_change=true, derivative=true)
 
     xlabel --> latexstring("t")
     labels = [LaTeXString("Net change"), LaTeXString("Time derivative")]
     fontfamily --> "Computer Modern"
 
-    @series begin
-        linestyle --> :solid
-        label --> labels[1]
-        results.t, results.E[:,e] .- first(results.E[:,e])
+    if net_change
+        @series begin
+            linestyle --> :solid
+            label --> labels[1]
+            results.t, results.E[:,e] .- first(results.E[:,e])
+        end
     end
-    @series begin
-        linestyle --> :dot
-        label --> labels[2]
-        results.t, results.dEdt[:,e]
+    if derivative
+        @series begin
+            linestyle --> :dot
+            label --> labels[2]
+            results.t, results.dEdt[:,e]
+        end
     end
 end

src/ConservationLaws/ConservationLaws.jl

@@ -123,7 +123,8 @@ module ConservationLaws
             initial_data::AbstractGridFunction{d};
             periodic::Bool=false) where {d, PDEType, N_c}
 
-            return new{d,typeof(conservation_law),typeof(initial_data),typeof(conservation_law.source_term)}(
+            return new{d,typeof(conservation_law),typeof(initial_data),
+                typeof(conservation_law.source_term)}(
                 conservation_law, initial_data, periodic, N_c)
         end
     end

src/ConservationLaws/euler_navierstokes.jl

@@ -220,50 +220,6 @@ end
     return SVector{d+2}(f_ρ, (f_ρ*V_avg[m] + p_avg*n[m] for m in 1:d)...,
         f_ρ*C + 0.5*(p_L*Vn_R + p_R*Vn_L))
 end
-
-
-"""
-Isentropic vortex problem, taken verbatim from the Trixi.jl examples (https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_euler_vortex.jl).
-
-Domain should be [-10,10] × [-10,10].
-"""
-struct TrixiIsentropicVortex <: AbstractGridFunction{2} 
-    γ::Float64
-    strength::Float64
-    N_c::Int
-    function TrixiIsentropicVortex(conservation_law::EulerEquations{2}, 
-        strength::Float64=5.0)
-        return new(conservation_law.γ,strength,4)
-    end
-end
-
-function evaluate(f::TrixiIsentropicVortex, x::NTuple{2,Float64},t::Float64=0.0)
-    inicenter = SVector(0.0, 0.0)
-    iniamplitude = f.strength
-
-    # base flow
-    gamma = f.γ
-    rho = 1.0
-    v1 = 1.0
-    v2 = 1.0
-    vel = SVector(v1, v2)
-    p = 25.0
-
-    rt = p / rho                  
-    t_loc = 0.0
-    cent = inicenter + vel*t_loc
-    cent = SVector(x) - cent 
-    cent = SVector(-cent[2], cent[1])
-    r2 = cent[1]^2 + cent[2]^2
-    du = iniamplitude / (2*π) * exp(0.5 * (1 - r2))
-    dtemp = -(gamma - 1) / (2 * gamma * rt) * du^2
-    rho = rho * (1 + dtemp)^(1 / (gamma - 1))
-    vel = vel + du * cent
-    v1, v2 = vel
-    p = p * (1 + dtemp)^(gamma / (gamma - 1))
-    return [rho, rho*v1, rho*v2, p/(gamma-1) + 0.5*rho*(v1^2 + v2^2)]
-end
-
 struct IsentropicVortex <: AbstractGridFunction{2} 
     γ::Float64
     Ma::Float64 
@@ -281,7 +237,8 @@ struct IsentropicVortex <: AbstractGridFunction{2}
     end
 end
 
-function evaluate(f::IsentropicVortex, x::NTuple{2,Float64},t::Float64=0.0)
+@inline function evaluate(f::IsentropicVortex, x::NTuple{2,Float64},
+    t::Float64=0.0)
     (; γ, Ma, θ, R, β, σ², x_0) = f
     x_rel = ((x[1] - x_0[1])/R, (x[2] - x_0[2])/R)
     Ω = β*exp(-0.5/σ² * (x_rel[1]^2 + x_rel[2]^2))
@@ -291,7 +248,7 @@ function evaluate(f::IsentropicVortex, x::NTuple{2,Float64},t::Float64=0.0)
     v = (Ma*cos(θ) + dv[1], Ma*sin(θ) + dv[2])
     p = (ρ^γ)/γ
     E = p/(γ-1) + 0.5*ρ*(v[1]^2 + v[2]^2)
-    return [ρ, ρ*v[1], ρ*v[2], E]
+    return SVector{4}(ρ, ρ*v[1], ρ*v[2], E)
 end
 
 """
@@ -314,11 +271,11 @@ struct EulerPeriodicTest{d} <: AbstractGridFunction{d}
     end
 end
 
-function evaluate(f::EulerPeriodicTest{d}, 
+@inline function evaluate(f::EulerPeriodicTest{d}, 
     x::NTuple{d,Float64},t::Float64=0.0) where {d}
 
     ρ = 1.0 + f.strength*sin(2π*sum(x[m] for m in 1:d)/f.L)
-    return [ρ, fill(ρ,d)..., 1.0/(f.γ-1.0) + 0.5*ρ*d]
+    return SVector{d+2}(ρ, fill(ρ,d)..., 1.0/(f.γ-1.0) + 0.5*ρ*d)
 end
 
 """Inviscid 3D Taylor-Green vortex, I think I got this version of it from Shadpey and Zingg, "Entropy-Stable Multidimensional Summation-by-Parts Discretizations on hp-Adaptive Curvilinear Grids for Hyperbolic Conservation Laws," JSC 2020.
@@ -335,13 +292,14 @@ struct TaylorGreenVortex <: AbstractGridFunction{3}
     end
 end
 
-function evaluate(f::TaylorGreenVortex,  x::NTuple{3,Float64}, t::Float64=0.0)
+@inline function evaluate(f::TaylorGreenVortex, 
+    x::NTuple{3,Float64}, t::Float64=0.0)
 
     p = 1.0/(f.γ * f.Ma^2) + 
         1.0/16 * (2*cos(2*x[1]) + 2*cos(2*x[2]) +
         cos(2*x[1])*cos(2*x[3]) + cos(2*x[2])*cos(2*x[3]))
     u = sin(x[1])*cos(x[2])*cos(x[3])
     v = -cos(x[1])*sin(x[2])*cos(x[3])
 
-    return [1.0, u, v, 0.0,  p/(f.γ-1.0) + 0.5*ρ*(u^2 + v^2)]
+    return SVector{5}(1.0, u, v, 0.0,  p/(f.γ-1.0) + 0.5*ρ*(u^2 + v^2))
 end

src/SpatialDiscretizations/SpatialDiscretizations.jl

@@ -13,7 +13,7 @@ module SpatialDiscretizations
     @reexport using StartUpDG: RefElemData, AbstractElemShape, Line, Quad, Tri, Tet, Hex, SBP
 
     export AbstractApproximationType, AbstractTensorProduct, AbstractMultidimensional,
-    NodalTensor, ModalTensor, ModalMulti, NodalMulti, ModalMultiDiagE, NodalMultiDiagE, AbstractReferenceMapping, AbstractMetrics, ExactMetrics, ChanWilcoxMetrics, NoMapping, ReferenceApproximation, GeometricFactors, SpatialDiscretization, apply_reference_mapping, reference_derivative_operators, check_normals, check_facet_nodes, check_sbp_property, centroids, trace_constant, dim, χ, warped_product
+    NodalTensor, ModalTensor, ModalMulti, NodalMulti, ModalMultiDiagE, NodalMultiDiagE, AbstractReferenceMapping, AbstractMetrics, ExactMetrics, ConservativeCurlMetrics, ChanWilcoxMetrics, NoMapping, ReferenceApproximation, GeometricFactors, SpatialDiscretization, apply_reference_mapping, reference_derivative_operators, check_normals, check_facet_nodes, check_sbp_property, centroids, trace_constant, dim, χ, warped_product
 
     abstract type AbstractApproximationType end
     abstract type AbstractTensorProduct <: AbstractApproximationType end
@@ -59,7 +59,8 @@ module SpatialDiscretizations
 
     abstract type AbstractMetrics end
     struct ExactMetrics <: AbstractMetrics end
-    struct ChanWilcoxMetrics <: AbstractMetrics end
+    struct ConservativeCurlMetrics <: AbstractMetrics end
+    const ChanWilcoxMetrics = ConservativeCurlMetrics
 
     """Operators for local approximation on reference element"""
     struct ReferenceApproximation{RefElemType, ApproxType, D_type, V_type,

src/SpatialDiscretizations/mesh.jl

@@ -15,8 +15,7 @@ struct ChanWarping{d} <: AbstractMeshWarping{d}
 end
 
 struct UniformWarping{d}  <: AbstractMeshWarping{d}
-    # Chan, Del Rey Fernandez, Carpenter (2019)
-    # Shadpey and Zingg (2020)
+    # Chan, Del Rey Fernandez, Carpenter (2019) or Shadpey and Zingg (2020)
     factor::Float64
     L::NTuple{d,Float64}
 end
@@ -276,10 +275,8 @@ function GeometricFactors(mesh::MeshData{d},
     return GeometricFactors(J_q, Λ_q, J_f, nJf, nJq)
 end
 
-
-function GeometricFactors(mesh::MeshData{2}, 
-    reference_element::RefElemData{2}, 
-    ::ChanWilcoxMetrics)
+function GeometricFactors(mesh::MeshData{2}, reference_element::RefElemData{2}, 
+    ::ConservativeCurlMetrics)
 
     (; x, y) = mesh
     (; nrstJ, Dr, Ds, Vq, Vf) = reference_element
@@ -337,8 +334,7 @@ function GeometricFactors(mesh::MeshData{2},
 end
 
 function GeometricFactors(mesh::MeshData{3}, 
-    reference_element::RefElemData{3}, 
-    ::ChanWilcoxMetrics)
+    reference_element::RefElemData{3,Hex}, ::ConservativeCurlMetrics)
 
     (; x, y, z) = mesh
     (; nrstJ, Dr, Ds, Dt, Vq, Vf) = reference_element
@@ -406,6 +402,96 @@ function GeometricFactors(mesh::MeshData{3},
     return GeometricFactors(J_q, Λ_q, J_f, nJf, nJq)
 end
 
+@views function GeometricFactors(mesh::MeshData{3}, 
+    reference_element::RefElemData{3,Tet}, ::ConservativeCurlMetrics)
+
+    (; x, y, z) = mesh
+    (; N, nrstJ) = reference_element
+
+    # need a degree N+1 element on which to compute the argument of the curl
+    # operator, which is approximated as a degree N+1 polynomial
+    reference_element_Nplus1 = RefElemData(Tet(),reference_element.N+1)
+    N_to_Nplus1 = vandermonde(Tet(), N, 
+        reference_element_Nplus1.r, reference_element_Nplus1.s,
+        reference_element_Nplus1.t) / reference_element.VDM
+    Nplus1_to_N = vandermonde(Tet(), N+1, 
+        reference_element.r, reference_element.s,
+        reference_element.t) / reference_element_Nplus1.VDM
+
+    # before evaluating metrics at quadrature nodes, they are brought 
+    # back to degree N interpolation nodes -- this is exact since they
+    # were lowered a degree by taking the curl
+    Vq = reference_element.Vq * Nplus1_to_N
+    Vf = reference_element.Vf * Nplus1_to_N
+
+    # note, here we assume that mesh is same N_q, N_f every element
+    N_q = size(mesh.xyzq[1],1)
+    N_f = size(mesh.xyzf[1],1)
+    N_e = size(mesh.xyzq[1],2)
+
+    # here we assume same number of nodes per face
+    N_fac = num_faces(reference_element.element_type)
+    nodes_per_face = N_f ÷ N_fac
+
+    J_q = Matrix{Float64}(undef, N_q, N_e)
+    J_f = Matrix{Float64}(undef, N_f, N_e)
+    nJf = Array{Float64, 3}(undef, 3, N_f, N_e)
+    nJq = Array{Float64, 4}(undef, 3, N_fac, N_q, N_e)
+    Λ_q = Array{Float64, 4}(undef, N_q, 3, 3, N_e)
+    Λ_f = Array{Float64, 4}(undef, N_f, 3, 3, N_e)
+
+    # Jacobian as degree N polynomial represented in degree N nodal basis
+    _, _, _, _, _, _, _, _, _, J = geometric_factors(x, y, z, 
+        reference_element.Dr, reference_element.Ds, reference_element.Dt)
+    mul!(J_q, reference_element.Vq, J)
+
+    # Metric terms as degree N polynomials represented in degree N+1 nodal basis
+    rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ, _ = geometric_factors(
+            N_to_Nplus1*x, N_to_Nplus1*y, N_to_Nplus1*z, 
+            reference_element_Nplus1.Dr, 
+            reference_element_Nplus1.Ds,
+            reference_element_Nplus1.Dt)
+
+    # Evaluate metric at volume quadrature nodes
+    mul!(Λ_q[:,1,1,:], Vq, rxJ)
+    mul!(Λ_q[:,2,1,:], Vq, sxJ)
+    mul!(Λ_q[:,3,1,:], Vq, txJ)
+    mul!(Λ_q[:,1,2,:], Vq, ryJ)
+    mul!(Λ_q[:,2,2,:], Vq, syJ)
+    mul!(Λ_q[:,3,2,:], Vq, tyJ)
+    mul!(Λ_q[:,1,3,:], Vq, rzJ)
+    mul!(Λ_q[:,2,3,:], Vq, szJ)
+    mul!(Λ_q[:,3,3,:], Vq, tzJ)
+
+    # Evaluate metric at volume facet quadrature nodes
+    mul!(Λ_f[:,1,1,:], Vf, rxJ)
+    mul!(Λ_f[:,2,1,:], Vf, sxJ)
+    mul!(Λ_f[:,3,1,:], Vf, txJ)
+    mul!(Λ_f[:,1,2,:], Vf, ryJ)
+    mul!(Λ_f[:,2,2,:], Vf, syJ)
+    mul!(Λ_f[:,3,2,:], Vf, tyJ)
+    mul!(Λ_f[:,1,3,:], Vf, rzJ)
+    mul!(Λ_f[:,2,3,:], Vf, szJ)
+    mul!(Λ_f[:,3,3,:], Vf, tzJ)
+
+    # compute normals
+    @inbounds for k in 1:N_e  
+        for i in 1:N_q, f in 1:N_fac
+            n_ref = Tuple(nrstJ[m][nodes_per_face*(f-1)+1] for m in 1:3)
+            for n in 1:3
+                nJq[n,f,i,k] = sum(Λ_q[i,m,n,k]*n_ref[m] for m in 1:3)
+            end
+        end
+        for i in 1:N_f
+            for m in 1:3
+                nJf[m,i,k] = sum(Λ_f[i,n,m,k]*nrstJ[n][i] for n in 1:3)
+            end
+            J_f[i,k] = sqrt(sum(nJf[m,i,k]^2 for m in 1:3))
+        end
+    end
+    return GeometricFactors(J_q, Λ_q, J_f, nJf, nJq)
+end
+
 function uniform_periodic_mesh(
     reference_approximation::ReferenceApproximation{<:RefElemData{3, Tet},
     <:Union{NodalTensor,ModalTensor}},

test/burgers_fluxdiff_1d.jl

@@ -7,6 +7,7 @@ function burgers_fluxdiff_1d()
 
     M = 20
     p = 7
+
     form = FluxDifferencingForm(
         inviscid_numerical_flux=EntropyConservativeNumericalFlux())
 
@@ -28,7 +29,7 @@ function burgers_fluxdiff_1d()
     h = L / (reference_approximation.N_p * spatial_discretization.N_e)
     dt = CFL * h / 1.0
 
-    sol = solve(ode_problem, CarpenterKennedy2N54(), adaptive=false, dt=dt,
+    solve(ode_problem, CarpenterKennedy2N54(), adaptive=false, dt=dt,
         save_everystep=false, callback=save_callback(results_path, (0.0,T),  
             floor(Int, T/(dt*50))))
 

test/euler_1d_gauss.jl

@@ -27,7 +27,7 @@ function euler_1d_gauss()
         inviscid_numerical_flux=EntropyConservativeNumericalFlux())
 
     ode = semidiscretize(conservation_law, spatial_discretization, exact_sol, 
-        form, (0.0, T), ReferenceOperator())
+        form, (0.0, T), ReferenceOperator(), BLASAlgorithm())
 
     results_path = save_project(conservation_law,
         spatial_discretization, exact_sol, form, (0.0, T),