Merge pull request #37 from ranocha/hr/weno

WENO5

.github/workflows/ci.yml

@@ -18,7 +18,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.5'
+          - '1.6'
           # - 'nightly'
         os:
           - ubuntu-latest

Project.toml

@@ -26,7 +26,7 @@ RecipesBase = "0.7, 0.8, 1.0"
 Reexport = "0.2, 1.0"
 Roots = "0.8, 1.0"
 StaticArrays = "1.0"
-julia = "1"
+julia = "1.6"
 
 [extras]
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"

src/HyperbolicDiffEq.jl

@@ -109,6 +109,7 @@ include("finite_volumes/numerical_fluxes.jl")
 include("finite_volumes/periodic_fv_1d.jl")
 include("finite_volumes/reconstructions.jl")
 include("finite_volumes/modified_eno.jl")
+include("finite_volumes/weno_jiang_shu.jl")
 include("finite_volumes/central_reconstruction.jl")
 
 include("tecno/tecno.jl")
@@ -162,6 +163,7 @@ export FirstOrderFV, UniformPeriodicReconstructedFV1D, ScalarUniformPeriodicTecn
 export CentralReconstruction
 export ModifiedENO, ClassicalChoiceENO, BiasedENOChoice,
        MinL2Choice, BiasedMinL2Choice, LexMinLegendreChoice
+export WENOJiangShu
 export GodunovFlux, LocalLaxFriedrichsFlux, HartenLaxVanLeerFlux, HLL, SuliciuFlux,
        KineticFlux
 export CentralFlux, MorinishiFlux, DucrosEtAlFlux, KennedyGruberFlux, PirozzoliFlux

src/finite_volumes/weno_jiang_shu.jl

@@ -0,0 +1,76 @@
+
+"""
+    WENOJiangShu{K}
+
+The classical weighted essentially non-oscillatory (WENO reconstruction of
+Jiang and Shu on `K` cells.
+"""
+struct WENOJiangShu{K, T<:Real} <: AbstractReconstruction
+  ε::T
+end
+
+WENOJiangShu{K}(ε) where {K} = WENOJiangShu{K, typeof(ε)}(ε)
+
+WENOJiangShu{K}() where {K} = WENOJiangShu{K}(1.0e-6)
+
+@inline stencil_width(::WENOJiangShu{K}) where {K} = K
+@inline stencil_width_val(::WENOJiangShu{K}) where {K} = Val{K}()
+
+@inline order(::WENOJiangShu{K}) where {K} = K
+
+
+function (weno::WENOJiangShu{5})(edge_u, cell, u_m2, u_m1, u_0, u_p1, u_p2, balance_law, meshx)
+    ε = weno.ε
+    T = typeof(ε)
+
+    # smoothness indicators of the three three-point stencils
+    β0 = 13 * (u_0  - 2 * u_p1 + u_p2)^2 / 12 + (3 * u_0 - 4 * u_p1 + u_p2)^2 / 4
+    β1 = 13 * (u_m1 - 2 * u_0  + u_p1)^2 / 12 + (u_m1 - u_p1)^2 / 4
+    β2 = 13 * (u_m2 - 2 * u_m1 + u_0 )^2 / 12 + (u_m2 - 4 * u_m1 + 3 * u_0)^2 / 4
+
+    # left edge
+    edge_u_l0 = ( 11*u_0 - 7*u_p1 + 2*u_p2 ) / 6
+    edge_u_l1 = ( 5*u_0 - u_p1 + 2*u_m1 ) / 6
+    edge_u_l2 = ( 2*u_0 + 5*u_m1 - u_m2 ) / 6
+
+    d_l0 = 1 * one(T) / 10
+    d_l1 = 3 * one(T) / 5
+    d_l2 = 3 * one(T) / 10
+
+    α_l0 = d_l0 / (ε + β0)^2
+    α_l1 = d_l1 / (ε + β1)^2
+    α_l2 = d_l2 / (ε + β2)^2
+
+    sum_α_l = α_l0 + α_l1 + α_l2
+
+    ω_l0 = α_l0 / sum_α_l
+    ω_l1 = α_l1 / sum_α_l
+    ω_l2 = α_l2 / sum_α_l
+
+    # right edge
+    edge_u_r0 = ( 2*u_0 + 5*u_p1 - u_p2 ) / 6
+    edge_u_r1 = ( 5*u_0 + 2*u_p1 - u_m1 ) / 6
+    edge_u_r2 = ( 11*u_0 - 7*u_m1 + 2*u_m2 ) / 6
+
+    d_r0 = 3 * one(T) / 10
+    d_r1 = 3 * one(T) / 5
+    d_r2 = 1 * one(T) / 10
+
+    α_r0 = d_r0 / (ε + β0)^2
+    α_r1 = d_r1 / (ε + β1)^2
+    α_r2 = d_r2 / (ε + β2)^2
+
+    sum_α_r = α_r0 + α_r1 + α_r2
+
+    ω_r0 = α_r0 / sum_α_r
+    ω_r1 = α_r1 / sum_α_r
+    ω_r2 = α_r2 / sum_α_r
+
+    # assign values
+    @inbounds begin
+        edge_u[1,cell] = ω_l0 * edge_u_l0 + ω_l1 * edge_u_l1 + ω_l2 * edge_u_l2
+        edge_u[2,cell] = ω_r0 * edge_u_r0 + ω_r1 * edge_u_r1 + ω_r2 * edge_u_r2
+    end
+
+    nothing
+end

test/modified_eno_test.jl

@@ -1,4 +1,6 @@
 using Test, OrdinaryDiffEq, HyperbolicDiffEq, Roots
+using Statistics: mean
+using LinearAlgebra: norm
 
 
 function calc_error(balance_law, uₐₙₐ, tmin, tmax, fnum, reconstruction, N, parallel)
@@ -90,6 +92,9 @@ println("BiasedMinL2Choice")
 @test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{4}(), BiasedMinL2Choice()), Ns_lo, Val{:serial}()) > 3.6
 @test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, ModifiedENO(Val{5}(), BiasedMinL2Choice()), Ns_lo, Val{:serial}()) > 4.8
 
+println("WENOJiangShu")
+@test calc_order_estimate(balance_law, uₐₙₐ, tspan, fnum, WENOJiangShu{5}(), Ns_lo, Val{:serial}()) > 4.8
+
 
 # Order reduction observed by Rogerson & Meiburg and Shu (1990).
 balance_law = ConstantLinearAdvection()