use cholesky for initial condition

Project.toml

@@ -18,7 +18,6 @@ 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/ConservationLaws/ConservationLaws.jl

@@ -45,16 +45,16 @@ module ConservationLaws
         @inbounds for i in axes(u_in, 1)
             f_s = compute_two_point_flux(conservation_law, two_point_flux,
                 u_in[i,:], u_out[i,:])
-            a = numerical_flux.λ*wave_speed(
+            a = 0.5*numerical_flux.λ*wave_speed(
                 conservation_law, u_in[i,:], u_out[i,:], 
                 Tuple(n[m][i] for m in 1:d))
             @inbounds for e in 1:N_c
                 temp = 0.0
                 du = u_out[i,e] - u_in[i,e]
                 @inbounds for m in 1:d
-                    temp = temp + f_s[e,m]*n[m][i]
+                    @muladd temp = temp + f_s[e,m]*n[m][i]
                 end
-                f_star[i,e] = temp - a * du
+                @muladd f_star[i,e] = temp - a * du
             end
         end
     end

src/GridFunctions/GridFunctions.jl

@@ -1,6 +1,6 @@
 module GridFunctions
 
-    export AbstractGridFunction, SumOfFunctions, ConstantFunction, InitialDataSine, InitialDataCosine, InitialDataGaussian, InitialDataGassner, BurgersSolution, SourceTermGassner, GaussianNoise, NoSourceTerm, evaluate
+    export AbstractGridFunction, SumOfFunctions, ConstantFunction, InitialDataSine, InitialDataCosine, InitialDataGaussian, InitialDataSinCos, DerivativeSinCos, InitialDataGassner, BurgersSolution, SourceTermGassner, GaussianNoise, NoSourceTerm, evaluate
 
     abstract type AbstractGridFunction{d} end
 
@@ -76,6 +76,29 @@ module GridFunctions
             return new(k,ϵ,1)
         end
     end
+
+
+    struct InitialDataSinCos <: AbstractGridFunction{2}
+        A::Float64  # amplitude
+        k::NTuple{2,Float64}  # wave number in each direction
+        N_c::Int
+
+        function InitialDataSinCos(A::Float64,
+            k::NTuple{2,Float64})
+            return new(A,k,1)
+        end
+    end
+
+    struct DerivativeSinCos <: AbstractGridFunction{2}
+        A::Float64  # amplitude
+        k::NTuple{2,Float64}  # wave number in each direction
+        N_c::Int
+
+        function DerivativeSinCos(A::Float64,
+            k::NTuple{2,Float64})
+            return new(A,k,1)
+        end
+end
     struct GaussianNoise{d} <: AbstractGridFunction{d}
         σ::Float64
         N_c::Int
@@ -142,7 +165,8 @@ module GridFunctions
         return [f.σ*randn() for e in 1:f.N_c]
     end
 
-    @inline function evaluate(f::AbstractGridFunction{d}, x::NTuple{d,Vector{Float64}}, t::Float64=0.0) where {d}
+    @inline function evaluate(f::AbstractGridFunction{d}, 
+        x::NTuple{d,AbstractVector{Float64}}, t::Float64=0.0) where {d}
         N = length(x[1])
         u0 = Matrix{Float64}(undef, N, f.N_c)
         @inbounds for i in 1:N
@@ -151,7 +175,9 @@ module GridFunctions
         return u0
     end
 
-    function evaluate(f::AbstractGridFunction{d}, x::NTuple{d,Matrix{Float64}},t::Float64=0.0) where {d}
+    function evaluate(f::AbstractGridFunction{d}, 
+        x::NTuple{d,AbstractMatrix{Float64}},
+        t::Float64=0.0) where {d}
         N, N_e = size(x[1])
         u0 = Array{Float64}(undef, N, f.N_c, N_e)
         @inbounds Threads.@threads for k in 1:N_e
@@ -160,7 +186,19 @@ module GridFunctions
         return u0
     end
 
-    function evaluate(::NoSourceTerm{d}, ::NTuple{d,Vector{Float64}}, ::Float64) where {d}
+    @inline function evaluate(::NoSourceTerm{d}, 
+        ::NTuple{d,Vector{Float64}}, ::Float64) where {d}
         return nothing
     end
-end
+
+    @inline function evaluate(f::InitialDataSinCos, 
+        x::NTuple{2,Float64},t::Float64=0.0)
+        return f.A*sin(f.k[1]*x[1])*cos(f.k[2]*x[2])
+    end
+
+    @inline function evaluate(f::DerivativeSinCos, 
+        x::NTuple{2,Float64},t::Float64=0.0)
+        return f.A*f.k[1]*cos(f.k[1]*x[1])*cos(f.k[2]*x[2])
+    end
+end
+

src/Solvers/preprocessing.jl

@@ -1,25 +1,25 @@
-function initialize(initial_data::AbstractGridFunction{d},
-    conservation_law::AbstractConservationLaw,
+@inline @views function initialize(initial_data::AbstractGridFunction{d},
     spatial_discretization::SpatialDiscretization{d}) where {d}
 
-    (; geometric_factors) = spatial_discretization
-    (; V, W) = spatial_discretization.reference_approximation
+    (; geometric_factors, N_e) = spatial_discretization
+    (; V, W, N_p) = spatial_discretization.reference_approximation
     (; xyzq) = spatial_discretization.mesh
-    N_p, N_c, N_e = get_dof(spatial_discretization, conservation_law)
+    (; N_c) = initial_data
+
+    u0 = Array{Float64}(undef, N_p, N_c, N_e)
 
     if V isa UniformScalingMap
-        u0 = Array{Float64}(undef, N_p, N_c, N_e)
         Threads.@threads for k in 1:N_e
-            u0[:,:,k] .= evaluate(initial_data, Tuple(xyzq[m][:,k] for m in 1:d),0.0)
+            u0[:,:,k] .= evaluate(
+                initial_data, Tuple(xyzq[m][:,k] for m in 1:d),0.0)
         end
     else
-        u0 = Array{Float64}(undef, N_p, N_c, N_e)
         Threads.@threads for k in 1:N_e
-            M =  Matrix(V' * W * Diagonal(geometric_factors.J_q[:,k]) * V)
-            rhs = similar(u0[:,:,k])
-            mul!(rhs, V' * W * Diagonal(geometric_factors.J_q[:,k]), 
+            M = Matrix(V' * W * Diagonal(geometric_factors.J_q[:,k]) * V)
+            factorization = cholesky(Symmetric(M))
+            mul!(u0[:,:,k], V' * W * Diagonal(geometric_factors.J_q[:,k]), 
                 evaluate(initial_data, Tuple(xyzq[m][:,k] for m in 1:d),0.0))
-            u0[:,:,k] .= M \ rhs
+            ldiv!(factorization, u0[:,:,k])
         end
     end
     return u0
@@ -46,7 +46,7 @@ function semidiscretize(
         end
     end
 
-    u0 = initialize(initial_data, conservation_law, spatial_discretization)
+    u0 = initialize(initial_data, spatial_discretization)
 
     if PDEType == SecondOrder && strategy isa ReferenceOperator
         @warn "Reference-operator approach only implemented for first-order equations. Using physical-operator formulation."