chan-wilcox metrics in 3D

examples/euler_vortex_2d.ipynb

@@ -1089,9 +1089,7 @@
     "error_analysis = ErrorAnalysis(results_path, conservation_law, \n",
     "    spatial_discretization)\n",
     "println(\"L2 error:\")\n",
-    "println(analyze(error_analysis, last(sol.u), exact_solution, T))\n",
-    "\n",
-    "#[0.02868402429755924, 0.04483188475445239, 0.041362524563684167, 0.05871228807076795]"
+    "println(analyze(error_analysis, last(sol.u), exact_solution, T))"
    ]
   },
   {

src/SpatialDiscretizations/mesh.jl

@@ -237,19 +237,14 @@ function GeometricFactors(mesh::MeshData{d},
 
     J_q = Matrix{Float64}(undef, N_q, N_e)
     J_f = Matrix{Float64}(undef, N_f, N_e)
-
     nJf = Array{Float64, 3}(undef, d, N_f, N_e)
     nJq = Array{Float64, 4}(undef, d, N_fac, N_q, N_e)
-
-    # first dimension is node index, 
-    # second and third are matrix indices mn,
-    # fourth is element index.
     dxdr_q = Array{Float64, 4}(undef, N_q, d, d, N_e)
     Λ_q = Array{Float64, 4}(undef, N_q, d, d, N_e)
     dxdr_f = Array{Float64, 4}(undef, N_f, d, d, N_e)   
 
-    for k in 1:N_e
-        @inbounds for m in 1:d, n in 1:d
+    @inbounds for k in 1:N_e
+        for m in 1:d, n in 1:d
             # evaluate metric at mapping nodes
             dxdr = reference_element.Drst[n]*mesh.xyz[m][:,k]
 
@@ -259,7 +254,7 @@ function GeometricFactors(mesh::MeshData{d},
         end
 
         # loops over slower indices
-        @inbounds for i in 1:N_q
+        for i in 1:N_q
             J_q[i,k], Λ_q[i,:,:,k] = metrics(SMatrix{d,d}(dxdr_q[i,:,:,k]))
             for f in 1:N_fac
                 n_ref = Tuple(nrstJ[m][nodes_per_face*(f-1)+1] for m in 1:d)
@@ -289,7 +284,6 @@ function GeometricFactors(mesh::MeshData{2},
     (; x, y) = mesh
     (; nrstJ, Dr, Ds, Vq, Vf) = reference_element
 
-
     # 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)
@@ -301,13 +295,8 @@ function GeometricFactors(mesh::MeshData{2},
 
     J_q = Matrix{Float64}(undef, N_q, N_e)
     J_f = Matrix{Float64}(undef, N_f, N_e)
-
     nJf = Array{Float64, 3}(undef, 2, N_f, N_e)
     nJq = Array{Float64, 4}(undef, 2, N_fac, N_q, N_e)
-
-    # first dimension is node index, 
-    # second and third are matrix indices mn,
-    # fourth is element index.
     Λ_q = Array{Float64, 4}(undef, N_q, 2, 2, N_e)
 
     @inbounds @views for k in 1:N_e
@@ -347,6 +336,76 @@ function GeometricFactors(mesh::MeshData{2},
     return GeometricFactors(J_q, Λ_q, J_f, nJf, nJq)
 end
 
+function GeometricFactors(mesh::MeshData{3}, 
+    reference_element::RefElemData{3}, 
+    ::ChanWilcoxMetrics)
+
+    (; x, y, z) = mesh
+    (; nrstJ, Dr, Ds, Dt, Vq, Vf) = reference_element
+
+    # 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)
+
+    @inbounds @views for k in 1:N_e
+
+        rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ, J = geometric_factors(
+            x[:,k], y[:,k], z[:,k], Dr, Ds, Dt)
+
+        mul!(Λ_q[:,1,1,k], Vq, rxJ)
+        mul!(Λ_q[:,2,1,k], Vq, sxJ)
+        mul!(Λ_q[:,3,1,k], Vq, txJ)
+        mul!(Λ_q[:,1,2,k], Vq, ryJ)
+        mul!(Λ_q[:,2,2,k], Vq, syJ)
+        mul!(Λ_q[:,3,2,k], Vq, tyJ)
+        mul!(Λ_q[:,1,3,k], Vq, rzJ)
+        mul!(Λ_q[:,2,3,k], Vq, szJ)
+        mul!(Λ_q[:,3,3,k], Vq, tzJ)
+        mul!(J_q[:,k], Vq, J)
+
+        Λ_f = Array{Float64, 3}(undef, N_f, 3, 3)
+        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)
+
+        # loops over slower indices
+        @inbounds for i in 1:N_q
+            for 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
+        end
+
+        # get scaled normal vectors - this includes scaling for ref. quadrature weights on long side of right-angled triangle.
+        @inbounds for i in 1:N_f
+            for m in 1:3
+                nJf[m,i,k] = sum(Λ_f[i,n,m]*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/advection_3d.jl

@@ -0,0 +1,55 @@
+function advection_3d()
+
+    a = (1.0,1.0,1.0)  # advection velocity
+    L = 1.0  # domain length
+    T = 1.0  # end time
+
+    conservation_law = LinearAdvectionEquation(a)
+    exact_solution = InitialDataCosine(1.0,(2π/L, 2π/L, 2π/L));
+
+    M = 2
+    p = 4
+
+    reference_approximation = ReferenceApproximation(
+        ModalTensor(p), Tet(), mapping_degree=4,
+        sum_factorize_vandermonde=false)
+
+    form = StandardForm(mapping_form=SkewSymmetricMapping(), 
+        inviscid_numerical_flux=CentralNumericalFlux())
+
+    uniform_mesh = uniform_periodic_mesh(reference_approximation,
+        ((0.0,L),(0.0,L),(0.0,L)), (M,M,M))
+
+    mesh = warp_mesh(uniform_mesh, reference_approximation, 0.1, L) 
+
+    spatial_discretization = SpatialDiscretization(mesh, 
+        reference_approximation, ChanWilcoxMetrics())
+
+    results_path = save_project(conservation_law,
+         spatial_discretization, exact_solution, form, (0.0, T),
+         "results/advection_3d/", overwrite=true, clear=true)
+    CFL = 0.1
+    h = L/(reference_approximation.N_p * spatial_discretization.N_e)^(1/3)
+    dt = CFL * h / sqrt(a[1]^2 + a[2]^2 + a[3]^2)
+
+    ode_problem = semidiscretize(conservation_law, spatial_discretization,
+        exact_solution, form, (0.0, T), ReferenceOperator(), BLASAlgorithm());
+
+    sol = solve(ode_problem, CarpenterKennedy2N54(), adaptive=false, dt=dt,
+        save_everystep=false, callback=save_callback(results_path, (0.0,T),
+        floor(Int, T/(dt*50))))
+
+    error_analysis = ErrorAnalysis(results_path, conservation_law, 
+        spatial_discretization, JaskowiecSukumarQuadrature(2p+3))
+
+    error_results = analyze(error_analysis, last(sol.u), exact_solution, T)
+
+    time_steps = load_time_steps(results_path)
+    conservation_results = analyze(PrimaryConservationAnalysis(results_path, 
+        conservation_law, spatial_discretization), time_steps)
+    energy_results = analyze(EnergyConservationAnalysis(results_path, 
+        conservation_law, spatial_discretization), time_steps)
+
+    return error_results[1], conservation_results.E[end,1] - conservation_results.E[1,1], maximum(abs.(energy_results.dEdt[:,1]))
+
+end

test/euler_vortex_2d_modal.jl

@@ -14,7 +14,7 @@ function euler_vortex_2d_modal(M::Int=4)
     p = 3
 
     form = FluxDifferencingForm(
-         inviscid_numerical_flux=EntropyConservativeNumericalFlux())
+         inviscid_numerical_flux=LaxFriedrichsNumericalFlux())
 
     reference_approximation = ReferenceApproximation(ModalTensor(p), 
         Tri(), mapping_degree=p)
@@ -44,14 +44,10 @@ function euler_vortex_2d_modal(M::Int=4)
     error_analysis = ErrorAnalysis(results_path, conservation_law, 
         spatial_discretization)
     error_results = analyze(error_analysis, last(sol.u), exact_solution, T)
-
     conservation_analysis = PrimaryConservationAnalysis(results_path, 
         conservation_law, spatial_discretization)
     conservation_results = analyze(conservation_analysis, 
         load_time_steps(results_path))
-    entropy_analysis = EntropyConservationAnalysis(results_path, 
-        conservation_law, spatial_discretization, ode.p.mass_solver)
-    entropy_results = analyze(entropy_analysis, load_time_steps(results_path))
-
-    return error_results, conservation_results.E[end,:] - conservation_results.E[1,:], maximum(abs.(entropy_results.dEdt[:,1]))
+
+    return error_results, conservation_results.E[end,:] - conservation_results.E[1,:]
 end

test/runtests.jl

@@ -6,9 +6,11 @@ include("burgers_fluxdiff_1d.jl")
 include("euler_1d_gauss.jl")
 include("euler_vortex_2d_nodal_diage.jl")
 include("euler_vortex_2d_modal.jl")
+include("advection_3d.jl")
 
 const tol = 1.0e-10
 const p = 4
+
 @testset "Advection-Diffusion 1D ModalMulti" begin
     (l2, conservation, _) = test_driver(
         ReferenceApproximation(ModalMulti(p), Line()), 
@@ -79,26 +81,17 @@ end
     @test entropy ≈ 0.0 atol=tol
 end
 
-@testset "Isentropic Euler vortex FluxDiff ModalTensor Tri 2D" begin
-    (l2, conservation, entropy) = euler_vortex_2d_modal()
+@testset "Isentropic Euler vortex FluxDiff ModalTensor Tri 2D LF" begin
+    (l2, conservation) = euler_vortex_2d_modal()
 
-    @test l2 ≈ [0.03734599714012373, 0.075883163090583, 0.06103456983950774, 0.07065687750221998] atol=tol
+    @test l2 ≈ [0.015568197027072704, 0.040539693811761104, 0.04060141777050208, 0.043960971468832745] atol=tol
     @test conservation ≈ [0.0, 0.0, 0.0, 0.0] atol=tol
-    @test entropy ≈ 0.0 atol=tol
 end
 
 @testset "Advection 3D Energy-Conservative ModalTensor Tet" begin
-    (l2, conservation, energy) = test_driver(
-        ReferenceApproximation(ModalTensor(p), Tet(), 
-            mapping_degree=3, sum_factorize_vandermonde=false),
-        LinearAdvectionEquation((1.0,1.0,1.0)),
-        InitialDataSine(1.0,(2*π, 2*π, 2*π)),
-        StandardForm(mapping_form=SkewSymmetricMapping(),
-            inviscid_numerical_flux=LaxFriedrichsNumericalFlux(0.0)),
-        ReferenceOperator(), BLASAlgorithm(), 1.0, 2, 1.0, 1.0/50.0, 0.1,
-        "test_advection_3d_collapsed_econ")
-
-    @test l2 ≈ 0.18604185425549752 atol=tol
+    (l2, conservation, energy) = advection_3d()
+
+    @test l2 ≈ 0.18698123719889992 atol=tol
     @test conservation ≈ 0.0 atol=tol
     @test energy ≈ 0.0 atol=tol
 end