small changes, try with dop853

Project.toml

@@ -18,6 +18,7 @@ Jacobi = "83f21c0b-4282-5fbc-9e3f-f6da3d2e584c"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 NodesAndModes = "7aca2e03-f7e2-4192-9ec8-f4ca66d597fb"

src/MatrixFreeOperators/tensor_product_2d.jl

@@ -1,4 +1,5 @@
-struct TensorProductMap2D{A_type,B_type,σᵢ_type,σₒ_type} <: LinearMaps.LinearMap{Float64}
+struct TensorProductMap2D{A_type,B_type,σᵢ_type,
+    σₒ_type} <: LinearMaps.LinearMap{Float64}
     A::A_type
     B::B_type
     σᵢ::σᵢ_type

src/MatrixFreeOperators/tensor_product_3d.jl

@@ -1,9 +1,10 @@
-struct TensorProductMap3D{A_type,B_type,C_type} <: LinearMaps.LinearMap{Float64}
+struct TensorProductMap3D{A_type,B_type,C_type,σᵢ_type, 
+    σₒ_type} <: LinearMaps.LinearMap{Float64}
     A::A_type
     B::B_type
     C::C_type
-    σᵢ::AbstractArray{Int,3}
-    σₒ::AbstractArray{Int,3}
+    σᵢ::σᵢ_type
+    σₒ::σₒ_type
 end
 
 @inline Base.size(L::TensorProductMap3D) = (
@@ -34,7 +35,7 @@ function TensorProductMap3D(A, B, C)
     elseif C isa Union{LinearMaps.WrappedMap,OctavianMap}
         C = SMatrix{M3,N3,Float64}(C.lmap)
     end
-    return TensorProductMap3D(A,B,C, σᵢ, σₒ)
+    return TensorProductMap3D(A, B, C, σᵢ, σₒ)
 end
 
 """

src/MatrixFreeOperators/warped_product_2d.jl

@@ -1,16 +1,16 @@
 """
 Warped tensor-product operator (e.g. for Dubiner-type bases)
 """    
-struct WarpedTensorProductMap2D{A_type,B_type,σᵢ_type,σₒ_type} <: LinearMaps.LinearMap{Float64}
+struct WarpedTensorProductMap2D{A_type,B_type,
+    σᵢ_type,σₒ_type} <: LinearMaps.LinearMap{Float64}
     A::A_type
     B::B_type
     σᵢ::σᵢ_type 
     σₒ::σₒ_type
     N2::Vector{Int}
     size::NTuple{2,Int}
 
-    function WarpedTensorProductMap2D(
-        A::A_type,B::B_type, σᵢ::σᵢ_type,
+    function WarpedTensorProductMap2D(A::A_type,B::B_type, σᵢ::σᵢ_type,
         σₒ::σₒ_type) where {A_type,B_type,σᵢ_type,σₒ_type}
         N2 = [count(a -> a>0, σᵢ[β1,:]) for β1 in axes(σᵢ,1)]
         return new{A_type,B_type,σᵢ_type,σₒ_type}(A,B,σᵢ,σₒ,
@@ -69,7 +69,7 @@ Evaluate the matrix-vector product
     (; A, B, σᵢ, σₒ, N2) = L.lmap
 
     Z = MMatrix{size(σᵢ,1), size(σₒ,2),Float64}(undef) # stack allocate
-    
+
     @inbounds for β1 in axes(σᵢ,1), α2 in axes(σₒ,2)
         temp = 0.0
         @simd for α1 in axes(σₒ,1)

src/Solvers/Solvers.jl

@@ -49,8 +49,6 @@ module Solvers
             LaxFriedrichsNumericalFlux()
         viscous_numerical_flux::ViscousNumericalFlux = BR1()
         two_point_flux::TwoPointFlux = EntropyConservativeFlux()
-        facet_correction::Bool = false
-        entropy_projection::Bool = false
     end
 
     struct ReferenceOperators{D_type, Dt_type, V_type, Vt_type,
@@ -64,9 +62,9 @@ module Solvers
         W::Diagonal{Float64, Vector{Float64}}
         B::Diagonal{Float64, Vector{Float64}}
         halfWΛ::Array{Diagonal{Float64, Vector{Float64}},3} # d x d x N_e
-        halfN::Matrix{Diagonal{Float64, Vector{Float64}}}
-        BJf::Vector{Diagonal{Float64, Vector{Float64}}}
-        n_f::Array{Float64,3}
+        halfN::Matrix{Diagonal{Float64, Vector{Float64}}} # d x N_e
+        BJf::Vector{Diagonal{Float64, Vector{Float64}}} # N_e
+        n_f::Array{Float64,3} # d x N_f x N_e
     end
 
     struct PhysicalOperators{d, VOL_type, FAC_type, V_type, 
@@ -75,7 +73,7 @@ module Solvers
         FAC::Vector{FAC_type}
         V::V_type
         R::R_type
-        n_f::Array{Float64,3}
+        n_f::Array{Float64,3} # d x N_f x N_e
     end
 
     struct FluxDifferencingOperators{S_type,
@@ -90,11 +88,11 @@ module Solvers
         W::Diagonal{Float64, Vector{Float64}}
         B::Diagonal{Float64, Vector{Float64}}
         WJ::Vector{Diagonal{Float64, Vector{Float64}}}
-        Λ_q::Array{Float64,4}
-        BJf::Vector{Diagonal{Float64, Vector{Float64}}}
-        n_f::Array{Float64,3}
-        halfnJf::Array{Float64,3}
-        halfnJq::Array{Float64,4}
+        Λ_q::Array{Float64,4} # N_q x d x d x N_e
+        BJf::Vector{Diagonal{Float64, Vector{Float64}}} # N_e
+        n_f::Array{Float64,3} # d x N_f x N_e
+        halfnJf::Array{Float64,3} # d x N_f x N_e
+        halfnJq::Array{Float64,4} # d x N_q x num_faces x N_e
         nodes_per_face::Int64
     end
 
@@ -230,7 +228,7 @@ module Solvers
     function Solver(
         conservation_law::AbstractConservationLaw{d,SecondOrder,N_c},     
         spatial_discretization::SpatialDiscretization{d},
-        form::StandardForm, ::PhysicalOperator, alg::AbstractOperatorAlgorithm,
+        form::StandardForm, ::AbstractStrategy, alg::AbstractOperatorAlgorithm,
         mass_solver::AbstractMassMatrixSolver, 
         parallelism::AbstractParallelism) where {d, N_c}
 
@@ -271,6 +269,7 @@ module Solvers
 
         @inbounds @views for k in 1:N_e
             WJ = Diagonal(W .* J_q[:,k])
+            # this will throw if M is not SPD
             M = cholesky(Symmetric(VDM' * WJ * VDM))
             lmul!(WJ, u_q[:,:,k])
             mul!(u0[:,:,k], VDM', u_q[:,:,k])

src/Solvers/flux_differencing_form.jl

@@ -70,18 +70,19 @@ end
     end
 end
 
+# no-op for LGL/diagonal-E operators
 @inline function facet_correction!(
-    r_q::AbstractMatrix{Float64}, # N_q x N_c 
-    f_f::AbstractMatrix{Float64}, # N_f x N_c
-    CORR::AbstractMatrix{Float64}, # N_f x N_q
-    conservation_law::AbstractConservationLaw{d},
-    two_point_flux::AbstractTwoPointFlux,
-    halfnJf::AbstractMatrix{Float64},
-    halfnJq::AbstractArray{Float64,3},
-    u_q::AbstractMatrix{Float64},
-    u_f::AbstractMatrix{Float64},
-    nodes_per_face::Int,
-    ::Val{false}) where {d}
+    ::AbstractMatrix{Float64}, 
+    ::AbstractMatrix{Float64}, 
+    ::Nothing, # no facet correction operator
+    ::AbstractConservationLaw{d},
+    ::AbstractTwoPointFlux,
+    ::AbstractMatrix{Float64},
+    ::AbstractArray{Float64,3},
+    ::AbstractMatrix{Float64},
+    ::AbstractMatrix{Float64},
+    ::Int) where {d}
+
     return
 end
 
@@ -95,8 +96,7 @@ end
     halfnJq::AbstractArray{Float64,3},
     u_q::AbstractMatrix{Float64},
     u_f::AbstractMatrix{Float64},
-    nodes_per_face::Int,
-    ::Val{true}) where {d, N_c}
+    nodes_per_face::Int) where {d, N_c}
 
     @inbounds for i in axes(u_q,1)
         for j in axes(u_f,1)
@@ -129,8 +129,7 @@ end
     halfnJq::AbstractArray{Float64,3},
     u_q::AbstractMatrix{Float64},
     u_f::AbstractMatrix{Float64},
-    nodes_per_face::Int,
-    ::Val{true}) where {d, N_c}
+    nodes_per_face::Int) where {d, N_c}
 
     C_nz = nonzeros(C)
     row_index = rowvals(C)
@@ -162,27 +161,27 @@ end
     end
 end
 
-# specialized for no entropy projection 
-@inline @views function entropy_projection!(
+# specialize for LGL/Diag-E nodal operators (no entropy projection)
+@inline function entropy_projection!(
     ::AbstractMassMatrixSolver,
-    ::AbstractConservationLaw,
+    conservation_law::AbstractConservationLaw,
     u_q::AbstractMatrix, 
     u_f::AbstractMatrix, 
     ::AbstractMatrix,
     ::AbstractMatrix,
     ::AbstractMatrix,
-    V::LinearMap,
+    V::UniformScalingMap, # nodal
     ::LinearMap,
-    R::LinearMap, ::Diagonal,
-    u::AbstractMatrix, ::Int,
-    ::Val{false})
+    R::SelectionMap, # diag-E
+    ::Diagonal,
+    u::AbstractMatrix, ::Int)
+
     mul!(u_q, V, u)
     mul!(u_f, R, u_q)
     return
 end
 
 # specialized for nodal schemes (not necessarily diagonal-E)
-# this is really an "entropy extrapolation" and not "projection"
 @inline @views function entropy_projection!(
     ::AbstractMassMatrixSolver,
     conservation_law::AbstractConservationLaw,
@@ -191,24 +190,24 @@ end
     w_q::AbstractMatrix,
     w_f::AbstractMatrix,
     w::AbstractMatrix,
-    V::UniformScalingMap, 
+    V::UniformScalingMap, # nodal
     ::LinearMap,
-    R::LinearMap, ::Diagonal,
-    u::AbstractMatrix, ::Int,
-    ::Val{true})
+    R::LinearMap, # not just SelectionMap
+    ::Diagonal,
+    u::AbstractMatrix,
+    ::Int)
 
     mul!(u_q, V, u)
     @inbounds for i in axes(u, 1)
-        conservative_to_entropy!(w_q[i,:], conservation_law,u_q[i,:])
+        conservative_to_entropy!(w_q[i,:], conservation_law, u_q[i,:])
     end
     mul!(w_f, R, w_q)
     @inbounds for i in axes(u_f, 1)
-        entropy_to_conservative!(u_f[i,:], conservation_law,w_f[i,:])
+        entropy_to_conservative!(u_f[i,:], conservation_law, w_f[i,:])
     end
 end
 
-# general (i.e. suitable for modal) approach
-# uses a full entropy projection
+# most general (i.e. suitable for modal) approach
 @inline @views function entropy_projection!(
     mass_solver::AbstractMassMatrixSolver,
     conservation_law::AbstractConservationLaw,
@@ -217,11 +216,12 @@ end
     w_q::AbstractMatrix,
     w_f::AbstractMatrix,
     w::AbstractMatrix,
-    V::LinearMap,
+    V::LinearMap, # not just UniformScalingMap
     Vᵀ::LinearMap,
-    R::LinearMap, WJ::Diagonal,
-    u::AbstractMatrix, k::Int,
-    ::Val{true})
+    R::LinearMap, # not just SelectionMap
+    WJ::Diagonal,
+    u::AbstractMatrix,
+    k::Int)
 
     # evaluate entropy variables in terms of nodal conservative variables
     mul!(u_q, V, u)
@@ -247,26 +247,40 @@ end
     end
 end
 
+# for scalar equations, no entropy projection
+@inline @views function nodal_values!(u::AbstractArray{Float64,3},
+    solver::Solver{<:AbstractConservationLaw{<:Any,1},
+    <:FluxDifferencingOperators, <:AbstractMassMatrixSolver,<:FluxDifferencingForm}, k::Int)
+
+    (; u_q, u_f) = solver.preallocated_arrays
+    (; V, R) = solver.operators
+
+    mul!(u_q[:,:,k], V, u[:,:,k])
+    mul!(u_f[:,k,:], R, u_q[:,:,k])
+
+    return
+end
+
+# for systems, dispatch entropy projection on V and R
 @inline @views function nodal_values!(u::AbstractArray{Float64,3},
     solver::Solver{<:AbstractConservationLaw,<:FluxDifferencingOperators,
     <:AbstractMassMatrixSolver,<:FluxDifferencingForm}, k::Int)
 
-    (; conservation_law, form, preallocated_arrays, mass_solver) = solver
+    (; conservation_law, preallocated_arrays, mass_solver) = solver
     (; f_f, u_q, r_q, u_f, temp) = preallocated_arrays
-    (; entropy_projection) = form
     (; V, Vᵀ, R, WJ) = solver.operators
 
     entropy_projection!(mass_solver, conservation_law, u_q[:,:,k], 
         u_f[:,k,:], r_q[:,:,k], f_f[:,:,k], temp[:,:,k], 
-        V, Vᵀ, R, WJ[k], u[:,:,k], k, Val(entropy_projection))
+        V, Vᵀ, R, WJ[k], u[:,:,k], k)
 end
 
 @inline @views function time_derivative!(dudt::AbstractArray{Float64,3},
     solver::Solver{<:AbstractConservationLaw,<:FluxDifferencingOperators,
     <:AbstractMassMatrixSolver,<:FluxDifferencingForm}, k::Int)
 
      (; conservation_law, connectivity, form, mass_solver) = solver
-     (; inviscid_numerical_flux, two_point_flux, facet_correction) = form
+     (; inviscid_numerical_flux, two_point_flux) = form
      (; f_f, u_q, r_q, u_f, CI) = solver.preallocated_arrays
      (; S, C, Vᵀ, Rᵀ, Λ_q, BJf, halfnJq, halfnJf, n_f, 
         nodes_per_face) = solver.operators
@@ -282,10 +296,10 @@ end
     flux_difference!(r_q[:,:,k], S, conservation_law, 
         two_point_flux, Λ_q[:,:,:,k], u_q[:,:,k])
 
-    # apply facet correction term (for operators w/o boundary nodes)
+    # apply facet correction term (if C is not nothing)
     facet_correction!(r_q[:,:,k], f_f[:,:,k], C, conservation_law,
         two_point_flux, halfnJf[:,:,k], halfnJq[:,:,:,k], 
-        u_q[:,:,k], u_f[:,k,:], nodes_per_face, Val(facet_correction))
+        u_q[:,:,k], u_f[:,k,:], nodes_per_face)
 
     # apply facet operators
     mul!(u_q[:,:,k], Rᵀ, f_f[:,:,k])