2D first draft

Author: Iddingsite

examples/1D_linear.jl examples/1D/1D_linear.jlexamples/1D_linear.jl renamed to examples/1D/1D_linear.jl

@@ -46,16 +46,16 @@ function main(nx=400)
 
 
     u = copy(u0_vec)
-    weno = WENOScheme(u; boundary=(2, 2), stag=true)
+    weno = WENOScheme(u; boundary=(2, 2), stag=true, multithreading=false)
 
     # advection velocity
-    a = ones(nx+1) .* 1
+    a = (;x=ones(nx+1))
 
     # grid size
     Δx = x[2] - x[1]
     Δt = CFL*Δx^(5/3)
 
-    tmax = period * (Lx + Δx) / maximum(abs.(a))
+    tmax = period * (Lx + Δx) / maximum(abs.(a.x))
 
     t = 0
 
@@ -76,4 +76,65 @@ function main(nx=400)
     display(f)
 end
 
-main()
+# main()
+
+
+nx = 400
+
+x_min = -1.0
+x_max = 1.0
+Lx = x_max - x_min
+
+x = range(x_min, stop=x_max, length=nx)
+
+# Courant number
+CFL = 0.4
+period = 4
+
+# Parameters for Shu test
+z = -0.7
+δ = 0.005
+β = log(2)/(36*δ^2)
+a = 0.5
+α = 10
+
+# Functions
+G(x, β, z) = exp.(-β .* (x .- z).^2)
+F(x, α, a) = sqrt.(max.(1 .- α^2 .* (x .- a).^2, 0.0))
+
+# Grid x assumed defined
+u0_vec = zeros(length(x))
+
+# Gaussian-like smooth bump at x in [-0.8, -0.6]
+idx = (x .>= -0.8) .& (x .<= -0.6)
+u0_vec[idx] .= (1/6) .* (G(x[idx], β, z - δ) .+ 4 .* G(x[idx], β, z) .+ G(x[idx], β, z + δ))
+
+# Heaviside step at x in [-0.4, -0.2]
+idx = (x .>= -0.4) .& (x .<= -0.2)
+u0_vec[idx] .= 1.0
+
+# Piecewise linear ramp at x in [0, 0.2]
+# Triangular spike at x=0.1, base width 0.2
+idx = abs.(x .- 0.1) .<= 0.1
+u0_vec[idx] .= 1 .- 10 .* abs.(x[idx] .- 0.1)
+
+# Elliptic/smooth bell at x in [0.4, 0.6]
+idx = (x .>= 0.4) .& (x .<= 0.6)
+u0_vec[idx] .= (1/6) .* (F(x[idx], α, a - δ) .+ 4 .* F(x[idx], α, a) .+ F(x[idx], α, a + δ))
+
+
+u = copy(u0_vec)
+weno = WENOScheme(u; boundary=(2, 2), stag=true, multithreading=false)
+
+# advection velocity
+a = (;x=ones(nx+1))
+
+# grid size
+Δx = x[2] - x[1]
+Δt = CFL*Δx^(5/3)
+
+tmax = period * (Lx + Δx) / maximum(abs.(a.x))
+
+t = 0
+@btime WENO_step!($u, $a, $weno, $Δt, $Δx)
+# 9.583 μs (0 allocations: 0 bytes)

examples/1D_linear_chmy.jl examples/1D/1D_linear_chmy.jlexamples/1D_linear_chmy.jl renamed to examples/1D/1D_linear_chmy.jl

@@ -78,8 +78,11 @@ function main(backend=CPU(), nx=200)
         end
     end
 
+
     f = plot(x, u0_vec, label="Exact", linestyle=:dash)
-    scatter!(f, x, u, label="WENO5", title="1D Linear Advection after $period periods", xlabel="x", ylabel="u", legend=:topright)
+    scatter!(f, x, u, label="WENO5", title="1D linear advection after $period periods", xlabel="x", ylabel="u", legend=:topright, markersize=2)
+    # save fig
+
     display(f)
 end
 

examples/2D/2D_linear.jl

@@ -0,0 +1,69 @@
+using FiniteDiffWENO5
+using Plots
+
+function main(;nx=400, ny=400)
+
+    Lx = 1.0
+    Δx = Lx / nx
+    Δy = Lx / ny
+
+    x = range(0, stop=Lx, length=nx)
+
+    # Courant number
+    CFL = 0.7
+    period = 1
+
+    # Grid x assumed defined
+    x = range(0, length=nx, stop= Lx)
+    y = range(0, length=ny, stop= Lx)
+    grid = (x .* ones(ny)', ones(nx) .* y')
+
+    w = π
+    vx0 = w .* (grid[1] .- Lx/2)
+    vy0 = -w .* (grid[2] .- Lx/2)
+
+    v = (;x=vy0, y=vx0)
+
+    x0 = 1/4
+    c = 0.08
+
+    u0 = zeros(ny, nx)
+
+    for I in CartesianIndices((ny, nx))
+        u0[I] = sign(exp(-((grid[1][I] - x0)^2 + (grid[2][I]' - x0)^2) / c^2) - 0.5) * 0.5 + 0.5
+    end
+
+    u = copy(u0)
+    weno = WENOScheme(u; boundary=(2, 2, 2, 2), stag=false, multithreading=true)
+
+
+    # grid size
+    Δt = CFL*min(Δx, Δy)^(5/3)
+
+    tmax = period / (w/(2*π))
+
+    t = 0
+    counter = 0
+
+    while t < tmax
+        WENO_step!(u, v, weno, Δt, Δx, Δy)
+
+
+        t += Δt
+
+        if t + Δt > tmax
+            Δt = tmax - t
+        end
+
+        if counter % 100 == 0
+            heatmap(x, y, u, title="t = $(round(t, digits=2))", clims=(0,1))
+            gui()
+        end
+
+        counter += 1
+
+    end
+
+end
+
+main(nx=400, ny=400)

ext/ChmyExt1D.jl

@@ -105,7 +105,7 @@ module ChmyExt1D
                     min(v.x[i], 0) * fr.x[i]
                     ) * Δx_
         else
-            du[i] = @muladd max(v[i], 0) * (fl.x[i+1] - fl.x[i]) * Δx_ + min(v[i], 0) * (fr.x[i+1] - fr.x[i]) * Δx_
+            du[i] = @muladd max(v.x[i], 0) * (fl.x[i+1] - fl.x[i]) * Δx_ + min(v.x[i], 0) * (fr.x[i+1] - fr.x[i]) * Δx_
         end
     end
 
@@ -122,7 +122,7 @@ module ChmyExt1D
     - `Δx`: The spatial grid size.
     - `grid::StructuredGrid`: The computational grid.
     """
-    function WENO_step!(u::T_field, v::NamedTuple{names, <:Tuple{<:T_field}}, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, grid::StructuredGrid, arch) where T_field <: AbstractVector{<:Real} where names
+    function WENO_step!(u::T_field, v::NamedTuple{(:x,), <:Tuple{<:T_field}}, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, grid::StructuredGrid, arch) where T_field <: AbstractVector{<:Real}
 
         backend = get_backend(u)
 

src/1D/semi_discretisation_1D.jl

@@ -1,39 +1,12 @@
 
 
-@inline function left_index(i, d, nx, ::Val{0})
-    # Dirichlet (clamped to domain)
-    return clamp(i - d, 1, nx)
-end
-
-@inline function left_index(i, d, nx, ::Val{1})
-    # Neumann (mirror the boundary value)
-    return max(i - d, 1)
-end
-
-@inline function left_index(i, d, nx, ::Val{2})
-    # Periodic (wrap around)
-    return mod1(i - d, nx)
-end
-
-@inline function right_index(i, d, nx, ::Val{0})
-    return clamp(i + d, 1, nx)   # Dirichlet
-end
-
-@inline function right_index(i, d, nx, ::Val{1})
-    return min(i + d, nx)        # Neumann
-end
-
-@inline function right_index(i, d, nx, ::Val{2})
-    return mod1(i + d, nx)       # Periodic
-end
-
 function WENO_flux!(fl, fr, u, weno, nx)
-    @unpack boundary, χ, γ, ζ, ϵ = weno
+    @unpack boundary, χ, γ, ζ, ϵ, multithreading = weno
 
     bL = Val(boundary[1])
     bR = Val(boundary[2])
 
-    @inbounds for i in 1:nx+1
+    @inbounds @maybe_threads multithreading for i in axes(fl.x,1)
         iwww = left_index(i, 3, nx, bL)
         iww  = left_index(i, 2, nx, bL)
         iw   = left_index(i, 1, nx, bL)
@@ -56,21 +29,21 @@ end
 
 function semi_discretisation_weno5!(du::T, v, weno::WENOScheme, Δx_) where T <: AbstractArray{<:Real, 1}
 
-    @unpack fl, fr, stag  = weno
+    @unpack fl, fr, stag, multithreading  = weno
 
     # use staggered grid or not for the velocities
     if stag
-        @inbounds for i = axes(du, 1)
+        @inbounds @maybe_threads multithreading for i = axes(du, 1)
             du[i] = @muladd (
-                            max(v[i+1], 0) * fl.x[i+1] +
-                            min(v[i+1], 0) * fr.x[i+1] -
-                            max(v[i], 0) * fl.x[i] -
-                            min(v[i], 0) * fr.x[i]
+                            max(v.x[i+1], 0) * fl.x[i+1] +
+                            min(v.x[i+1], 0) * fr.x[i+1] -
+                            max(v.x[i], 0) * fl.x[i] -
+                            min(v.x[i], 0) * fr.x[i]
                             ) * Δx_
         end
     else
-        @inbounds for i = axes(du, 1)
-            du[i] = @muladd max(v[i], 0) * (fl.x[i+1] - fl.x[i]) * Δx_ + min(v[i], 0) * (fr.x[i+1] - fr.x[i]) * Δx_
+        @inbounds @maybe_threads multithreading for i = axes(du, 1)
+            du[i] = @muladd max(v.x[i], 0) * (fl.x[i+1] - fl.x[i]) * Δx_ + min(v.x[i], 0) * (fr.x[i+1] - fr.x[i]) * Δx_
         end
     end
 end

src/1D/time_stepping.jl

@@ -16,26 +16,26 @@ function WENO_step!(u::T, v, weno::WENOScheme, Δt, Δx) where T <: AbstractVect
     nx = size(u, 1)
     Δx_ = inv(Δx)
 
-    @unpack ut, du, stag, fl, fr = weno
+    @unpack ut, du, stag, fl, fr, multithreading = weno
 
     WENO_flux!(fl, fr, u, weno, nx)
     semi_discretisation_weno5!(du, v, weno, Δx_)
 
-    @inbounds for i = axes(ut, 1)
+    @inbounds @maybe_threads multithreading for i = axes(ut, 1)
         ut[i] = @muladd u[i] - Δt * du[i]
     end
 
     WENO_flux!(fl, fr, ut, weno, nx)
     semi_discretisation_weno5!(du, v, weno, Δx_)
 
-    @inbounds for i = axes(ut, 1)
+    @inbounds @maybe_threads multithreading for i = axes(ut, 1)
         ut[i] = @muladd 0.75 * u[i] + 0.25 * ut[i] - 0.25 * Δt * du[i]
     end
 
     WENO_flux!(fl, fr, ut, weno, nx)
     semi_discretisation_weno5!(du, v, weno, Δx_)
 
-    @inbounds for i = axes(u, 1)
+    @inbounds @maybe_threads multithreading for i = axes(u, 1)
         u[i] = @muladd 1.0/3.0 * u[i] + 2.0/3.0 * ut[i] - (2.0/3.0) * Δt * du[i]
     end