entropy conservation analysis

src/Analysis/Analysis.jl

@@ -1,7 +1,7 @@
 module Analysis
 
     using LinearMaps: LinearMap, UniformScalingMap
-    using LinearAlgebra: Diagonal, dot, norm, eigen, inv, svd, qr, pinv, eigsortby, I, rank, cond
+    using LinearAlgebra: Diagonal, mul!, lmul!, dot, norm, eigen, inv, svd, qr, pinv, eigsortby, I, rank, cond
     using JLD2: save, load, save_object, load_object
     using Plots: plot, savefig, plot!, scatter, text, annotate!, vline!, grid, theme_palette, twinx, @layout
     using RecipesBase
@@ -15,10 +15,10 @@ module Analysis
     using Markdown
     using TimerOutputs
 
-    using ..ConservationLaws: AbstractConservationLaw
+    using ..ConservationLaws: AbstractConservationLaw, entropy, conservative_to_entropy
     using ..SpatialDiscretizations: SpatialDiscretization, ReferenceApproximation, uniform_periodic_mesh, quadrature
     using ..GridFunctions: AbstractGridFunction, evaluate
-    using ..Solvers: AbstractResidualForm, AbstractStrategy, Solver, AbstractMassMatrixSolver, WeightAdjustedSolver, mass_matrix, semidiscretize, LinearResidual, get_dof, rhs!
+    using ..Solvers: AbstractResidualForm, AbstractStrategy, Solver, AbstractMassMatrixSolver, WeightAdjustedSolver, mass_matrix,  mass_matrix_inverse, mass_matrix_solve!, semidiscretize, LinearResidual, get_dof, rhs!
     using ..File: new_path, load_project, load_solution, load_time_steps, load_snapshots, load_snapshots_with_derivatives, load_solver, save_callback, save_solution, save_project
     using ..Visualize: Plotter
 
@@ -40,7 +40,7 @@ module Analysis
     export LinearAnalysis, DynamicalAnalysisResults, KoopmanAnalysis, AbstractKoopmanAlgorithm, StandardDMD, ExtendedDMD, KernelDMD, KernelResDMD, ExtendedResDMD, GeneratorDMD, AbstractSamplingAlgorithmx, GaussianSampling, analyze_running, forecast, monomial_basis, monomial_derivatives, make_dmd_matrices, dmd, generate_samples
     include("dynamics.jl")
 
-    export ConservationAnalysis, PrimaryConservationAnalysis, EnergyConservationAnalysis, ConservationAnalysisResults, ConservationAnalysisResultsWithDerivative, plot_evolution
+    export ConservationAnalysis, PrimaryConservationAnalysis, EnergyConservationAnalysis, EntropyConservationAnalysis, ConservationAnalysisResults, ConservationAnalysisResultsWithDerivative, plot_evolution
     include("conservation.jl")
 
     export RefinementAnalysis, RefinementErrorAnalysis, RefinementAnalysisResults, RefinementErrorAnalysisResults, run_refinement, get_tickslogscale

src/Analysis/conservation.jl

@@ -25,6 +25,17 @@ struct EnergyConservationAnalysis <: ConservationAnalysis
     dict_name::String
 end
 
+struct EntropyConservationAnalysis <: ConservationAnalysis
+    mass_solver::AbstractMassMatrixSolver
+    WJ::Vector{Diagonal}
+    conservation_law::AbstractConservationLaw
+    N_c::Int
+    N_e::Int
+    V::LinearMap
+    results_path::String
+    dict_name::String
+end
+
 struct ConservationAnalysisResults <:AbstractConservationAnalysisResults
     t::Vector{Float64}
     E::Matrix{Float64}
@@ -58,9 +69,8 @@ function EnergyConservationAnalysis(results_path::String,
 
     _, N_c, N_e = get_dof(spatial_discretization, conservation_law)
 
-    (; W, V) =  spatial_discretization.reference_approximation
-    (; mesh, N_e) = spatial_discretization
-    (; J_q) = spatial_discretization.geometric_factors
+    (; V) =  spatial_discretization.reference_approximation
+    (; N_e) = spatial_discretization
 
     return EnergyConservationAnalysis(
         mass_solver, N_c, N_e, V ,results_path, "energy.jld2")
@@ -81,6 +91,22 @@ function EnergyConservationAnalysis(
         spatial_discretization, mass_solver)
 end
 
+function EntropyConservationAnalysis(results_path::String,
+    conservation_law::AbstractConservationLaw, 
+    spatial_discretization::SpatialDiscretization{d},
+    mass_solver=WeightAdjustedSolver(spatial_discretization)) where {d}
+
+    _, N_c, N_e = get_dof(spatial_discretization, conservation_law)
+
+    (; W, V) =  spatial_discretization.reference_approximation
+    (; geometric_factors, mesh, N_e) = spatial_discretization
+
+    WJ = [Diagonal(W .* geometric_factors.J_q[:,k]) for k in 1:N_e]
+
+    return EntropyConservationAnalysis(mass_solver, WJ, conservation_law, 
+        N_c, N_e, V, results_path, "entropy.jld2")
+end
+
 function evaluate_conservation(
     analysis::PrimaryConservationAnalysis, 
     u::Array{Float64,3})
@@ -90,6 +116,20 @@ function evaluate_conservation(
         for k in 1:N_e) for e in 1:N_c]
 end
 
+function evaluate_conservation(
+    analysis::EntropyConservationAnalysis, 
+    u::Array{Float64,3})
+    (; WJ, conservation_law, N_e, V) = analysis 
+    S = 0.0
+    @views for k in 1:N_e
+        u_q = V*u[:,:,k]
+        for i in axes(u,1)
+            S += WJ[k][i,i]*entropy(conservation_law, u_q[i,:])
+        end
+    end
+    return [S]
+end
+
 function evaluate_conservation(
     analysis::EnergyConservationAnalysis, 
     u::Array{Float64,3})
@@ -132,11 +172,62 @@ function evaluate_conservation_residual(
     return dEdt
 end
 
+function evaluate_conservation_residual(
+    analysis::EntropyConservationAnalysis, 
+    u::Array{Float64,3},
+    dudt::Array{Float64,3})
+    (; mass_solver, conservation_law, N_c, N_e, V, WJ) = analysis 
+
+    u_q = Matrix{Float64}(undef,size(V,1),N_c)
+    w_q = Matrix{Float64}(undef,size(V,1),N_c)
+    w = Matrix{Float64}(undef,size(u,1),N_c)
+    r = Matrix{Float64}(undef,size(u,1),N_c)
+    dEdt = 0.0
+    for k in 1:N_e
+        M = mass_matrix(mass_solver, k)
+        mul!(u_q,V,u[:,:,k])
+        for i in axes(u,1)
+            w_q[i,:] .= conservative_to_entropy(conservation_law, u_q[i,:])
+        end
+        P = mass_matrix_inverse(mass_solver, k) * V' * WJ[k]
+        for e in 1:N_c
+            dEdt += (P*w_q[:,e])'*M*dudt[:,e,k]
+        end
+    end
+    return [dEdt]
+
+end
+
+function analyze(analysis::EntropyConservationAnalysis,
+    time_steps::Vector{Int})
+
+    (; results_path, dict_name) = analysis
+    N_t = length(time_steps)
+    t = Vector{Float64}(undef,N_t)
+    E = Matrix{Float64}(undef,N_t, 1)
+    dEdt = Matrix{Float64}(undef,N_t, 1)
+
+    for i in 1:N_t
+        u, dudt, t[i] = load_solution(results_path, time_steps[i], load_du=true)
+        E[i,:] = evaluate_conservation(analysis, u)
+        dEdt[i,:] = evaluate_conservation_residual(analysis, u, dudt)
+    end
+
+    results = ConservationAnalysisResults(t,E)
+
+    save(string(results_path, dict_name), 
+    Dict("conservation_analysis" => analysis,
+        "conservation_results" => results))
+
+    return ConservationAnalysisResultsWithDerivative(t, E, dEdt)
+end
+
+
 function analyze(analysis::ConservationAnalysis, 
     initial_time_step::Union{Int,String}=0, 
     final_time_step::Union{Int,String}="final")
 
-    (; results_path, dict_name) = analysis
+    (; results_path) = analysis
 
     u_0, _ = load_solution(results_path, initial_time_step)
     u_f, _ = load_solution(results_path, final_time_step)
@@ -246,7 +337,6 @@ end
         label --> labels[2]
         results.t, results.dEdt[:,e]
     end
-
 end
 
 function plot_evolution(analysis::ConservationAnalysis, 

src/ConservationLaws/ConservationLaws.jl

@@ -5,7 +5,7 @@ module ConservationLaws
     using LinearAlgebra: mul!, I
     import ..GridFunctions: AbstractGridFunction, NoSourceTerm, InitialDataSine, InitialDataGaussian, InitialDataGassner, SourceTermGassner, evaluate
 
-    export physical_flux!, numerical_flux!, conservative_to_primitive, conservative_to_entropy, entropy_to_conservative, compute_two_point_flux, wave_speed, logmean, AbstractConservationLaw, AbstractPDEType, FirstOrder, SecondOrder, AbstractInviscidNumericalFlux, AbstractViscousNumericalFlux, NoInviscidFlux, NoViscousFlux, LaxFriedrichsNumericalFlux, RoeNumericalFlux, BR1, EntropyConservativeNumericalFlux, AbstractTwoPointFlux, EntropyConservativeFlux, NoTwoPointFlux, ExactSolution
+    export physical_flux!, numerical_flux!, entropy, conservative_to_primitive, conservative_to_entropy, entropy_to_conservative, compute_two_point_flux, wave_speed, logmean, AbstractConservationLaw, AbstractPDEType, FirstOrder, SecondOrder, AbstractInviscidNumericalFlux, AbstractViscousNumericalFlux, NoInviscidFlux, NoViscousFlux, LaxFriedrichsNumericalFlux, RoeNumericalFlux, BR1, EntropyConservativeNumericalFlux, AbstractTwoPointFlux, EntropyConservativeFlux, NoTwoPointFlux, ExactSolution
 
     abstract type AbstractConservationLaw{d, PDEType} end
     abstract type AbstractPDEType end
@@ -83,12 +83,11 @@ module ConservationLaws
         Algorithm based on the Taylor series trick from Ismail and Roe (2009). There are further optimizations that could be made, but let's leave it like this for now.
     """
     @inline function logmean(x::Float64, y::Float64)
-        # f = (y/x - 1) / (y/x + 1) = (y - x) / (x + y)
-        # rearrange to avoid divisions 
-        # uses trick from https://trixi-framework.github.io/Trixi.jl/stable/reference-trixi/#Trixi.ln_mean
+        # f = (y/x - 1) / (y/x + 1)
+        #   = (y - x) / (x + y)
+        # rearrange to avoid divisions using trick from https://trixi-framework.github.io/Trixi.jl/stable/reference-trixi/#Trixi.ln_mean
         f² = (x * (x - 2 * y) + y * y) / (x * (x + 2 * y) + y * y)
         if f² < 1.0e-4
-            println("using Ismail-roe trick")
             return (x + y) * 105 / (210 + f² * (70 + f² * (42 + f² * 30)))
             # faster factorized way to compute
             # (x + y) / (2 + 2/3 * f^2 + 2/5 * f^4 + 2/7 * f^6)

src/ConservationLaws/burgers.jl

@@ -104,6 +104,8 @@ function numerical_flux!(f_star::AbstractMatrix{Float64},
 end
 
 @inline conservative_to_primitive(::BurgersType, u) = u
+@inline entropy(::BurgersType, u) = 0.5*u[1]^2
+@inline conservative_to_entropy(::BurgersType, u) = u
 
 @inline function wave_speed(conservation_law::BurgersType{d},
     u_in::AbstractVector{Float64}, u_out::AbstractVector{Float64},

src/ConservationLaws/euler_navierstokes.jl

@@ -48,13 +48,6 @@ end
 
 const EulerType{d} = Union{EulerEquations{d}, NavierStokesEquations{d}}
 
-@inline function conservative_to_primitive(conservation_law::EulerType{d},
-    u::AbstractVector{Float64}) where {d}
-    return vcat(u[1], SVector{d}(u[m+1] / u[1] for m in 1:d),
-        (conservation_law.γ-1) * (u[end] - (0.5/u[1]) * 
-        (sum(u[m+1]^2 for m in 1:d))))
-end
-
 """
 Evaluate the flux for the Euler equations
 """
@@ -73,6 +66,29 @@ function physical_flux!(f::AbstractArray{Float64,3},
     end
 end
 
+@inline function entropy(conservation_law::EulerType{d}, 
+    u::AbstractVector{Float64}) where {d}
+    p = (conservation_law.γ-1) * (u[end] - (0.5/u[1]) * 
+        (sum(u[m+1]^2 for m in 1:d)))
+    return -u[1]*log(p/u[1]^conservation_law.γ)/(conservation_law.γ-1)
+end
+
+@inline function conservative_to_primitive(conservation_law::EulerType{d},
+    u::AbstractVector{Float64}) where {d}
+    return vcat(u[1], SVector{d}(u[m+1] / u[1] for m in 1:d),
+        (conservation_law.γ-1) * (u[end] - (0.5/u[1]) * 
+        (sum(u[m+1]^2 for m in 1:d))))
+end
+
+@inline function conservative_to_entropy(conservation_law::EulerType{d}, 
+    u::AbstractVector{Float64}) where {d}
+    (; γ) = conservation_law
+    k = (0.5/u[1]) * (sum(u[m+1]^2 for m in 1:d))
+    p = (γ-1) * (u[end] - k)
+    return vcat((γ-log(p/(u[1]^γ)))/(γ-1) - k/p,
+        SVector{d}(u[m+1] / p for m in 1:d), -u[1]/p)
+end
+
 @inline function wave_speed(conservation_law::EulerType{d},
     u_in::AbstractVector{Float64}, u_out::AbstractVector{Float64},
     n_f::NTuple{d, Float64}) where {d}

src/ConservationLaws/linear_advection_diffusion.jl

@@ -111,6 +111,8 @@ function numerical_flux!(f_star::AbstractMatrix{Float64},
 end
 
 @inline conservative_to_primitive(::AdvectionType, u) = u
+@inline conservative_to_entropy(::AdvectionType, u) = u
+@inline entropy(::AdvectionType, u) = 0.5*u[1]^2
 
 @inline function wave_speed(conservation_law::AdvectionType{d},
     ::AbstractVector{Float64}, ::AbstractVector{Float64},