runic formatting

.github/workflows/Runic.yml

@@ -12,7 +12,7 @@ jobs:
     name: Runic
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
       - uses: fredrikekre/runic-action@v1
         with:
           version: '1'

README.md

@@ -99,4 +99,4 @@ end
 
 Which produces the following result:
 
-![](/docs/assets/1D_linear_advection.png)
+![](/docs/src/assets/1D_linear_advection.png)

docs/make.jl

@@ -1,5 +1,15 @@
-push!(LOAD_PATH,"../src/")
+push!(LOAD_PATH, "../src/")
 
 using Documenter, FiniteDiffWENO5
 
-makedocs(sitename="My Documentation", remotes = nothing)
+makedocs(
+    sitename = "FiniteDiffWENO5",
+    remotes = nothing,
+    authors = "Hugo Dominguez",
+    pages = [
+        "Home" => "index.md",
+        "Getting Started" => "GettingStarted.md",
+        "Background" => "background.md",
+        "API" => "API.md",
+    ],
+)

docs/src/API.md

@@ -0,0 +1,6 @@
+# [API](@id API)
+
+```@autodocs
+Modules  = [FiniteDiffWENO5]
+
+```

docs/src/GettingStarted.md

@@ -0,0 +1,99 @@
+
+## Getting Started
+
+There are currently two exported functions in FiniteDiffWENO5.jl: `WENOScheme()` and `WENO_step!()`. The user must define the grid and the initial conditions themselves.
+
+`WENOScheme()` is used to create a WENO scheme structure containing all the necessary information for the WENO method, while `WENO_step!()` performs one step of the time integration using the WENO-Z method and a 3rd-order Runge-Kutta method. Refer to the docstrings to see the available options for each function.
+
+Here is a simple example of how to use FiniteDiffWENO5.jl to solve the 1D advection equation using the conservative form on a staggered grid. For more examples, please refer to the folder examples in the repository or the tests in the test folder.
+
+```julia
+using FiniteDiffWENO5
+
+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)
+
+# fix here boundary to periodic when equal to 2
+# staggering = true means that the advection velocity is defined on the sides of the cells and should be of size nx+1 compared to the scalar field u.
+# Multithreading to false means that the computations will be done on a single thread.
+# Set to true to enable multithreading if the Julia session was started with multiple threads (e.g. `julia -t 4`).
+weno = WENOScheme(u; boundary = (2, 2), stag = true, multithreading = false)
+
+# advection velocity with size nx+1 for staggered grid (here constant)
+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
+
+while t < tmax
+    # take as input the scalar field u, the advection velocity a as a NamedTuple,
+    # the WENO scheme struct weno, the time step Δt and the grid size Δx
+    WENO_step!(u, a, weno, Δt, Δx)
+
+    t += Δt
+
+    if t + Δt > tmax
+        Δt = tmax - t
+    end
+end
+
+f = Figure(size = (800, 600), dpi = 400)
+ax = Axis(f[1, 1], title = "1D linear advection after $period periods", xlabel = "x", ylabel = "u")
+lines!(ax, x, u0_vec, label = "Exact", linestyle = :dash, color = :red)
+scatter!(ax, x, u, label = "WENO5")
+xlims!(ax, x_min, x_max)
+axislegend(ax)
+display(f)
+```
+
+which outputs:
+
+![1D advection](./assets/1D_linear_advection.png)

docs/src/background.md

@@ -0,0 +1,20 @@
+# [Background](@id background)
+
+The weighted essentially non-oscillatory (WENO) schemes form a class of high-order accurate numerical methods for solving hyperbolic partial differential equations (PDEs). They are particularly effective at resolving sharp gradients and discontinuities without introducing spurious oscillations. The WENO schemes achieve this by constructing nonlinear convex combinations of lower-order candidate polynomials, where the weights are determined by the local smoothness of the solution. This approach allows the method to retain high-order accuracy in smooth regions while automatically reducing to lower-order, more stable behavior near discontinuities. The WENO framework was first introduced by [Liu et al. 1994](https://doi.org/10.1006/jcph.1994.1187]), building upon the earlier essentially non-oscillatory (ENO) schemes developed by [Harten et al. 1987](https://doi.org/10.1016/0021-9991(87)90031-3).
+
+In practice, WENO schemes can be formulated either in a finite-volume or finite-difference framework. In this package, we adopt the finite-difference formulation, which is particularly well-suited for problems defined on structured grids due to its simplicity and computational efficiency. The time integration is based on a third order strong stability preserving Runge-Kutta (SSP-RK3) method.
+
+The implementation of a finite-difference WENO scheme involves the following main steps:
+
+2. **Smoothness Indicators**: For each candidate stencil, a smoothness indicator is computed. This indicator quantifies how smooth the approximation is within that stencil, with lower values indicating smoother regions.
+3. **Weight Calculation**: Nonlinear weights are computed based on the smoothness indicators. Stencils with lower smoothness indicators receive higher weights, allowing the scheme to adaptively favor smoother regions.
+4. **Reconstruction**: The final high-order approximation is obtained by combining the candidate polynomials using the computed weights.
+5. **Flux Evaluation**: The reconstructed values are used to approximate the variable of interest at the cell interfaces, which are then used in the numerical flux calculations.
+
+In FiniteDiffWENO5.jl, two forms of advection operators are currently supported:
+1. **Non-conservative Form** $\mathbf{v} \cdot \nabla u$, where the velocity field $\mathbf{v}$ and scalar field $u$ are both defined at the same grid locations (collocated grid), and
+2. **Conservative Form** $\nabla \cdot (\mathbf{v}u)$, where $\mathbf{v}$ is defined on cell faces and $u$ at cell centers (staggered grid).
+
+In both formulations, only the scalar field $u$ is reconstructed at the cell interfaces using the WENO scheme.
+
+The package currently implements the WENO-Z reconstruction developed by [Borges et al. (2008)](https://doi.org/10.1016/j.jcp.2007.11.038). This variant introduces a modified computation of the nonlinear weights that improves accuracy near critical points—where the first derivative of the solution vanishes—while preserving the robust, non-oscillatory behavior of the classical WENO methods. Additional reconstruction variants may be included in future versions.

docs/src/index.md

@@ -1,8 +1,17 @@
+# [FiniteDiffWENO5.jl](@id home)
 
-```@index
-```
+FiniteDiffWENO5.jl is a Julia package that implements fifth-order finite-difference weighted essentially non-oscillatory (WENO) schemes for solving hyperbolic partial differential equations (PDEs) in 1D, 2D and 3D on regular grids.
+
+The package currently focuses on the non-conservative form of the advection terms ($\mathbf{v} \cdot \nabla u$) on a collocated grid, and the conservative form ($\nabla \cdot$ ($\mathbf{v} u$)) where the velocity field $\mathbf{v}$ and scalar field $u$ are on a staggered grid with the advection velocity located on the sides of the cells.
+
+The core of the package is written in pure Julia, focusing on performance using CPUs, but GPU support is available using KernelAbstractions.jl and Chmy.jl via an extension.
 
-```@autodocs
-Modules  = [FiniteDiffWENO5]
+## Installation
 
-```
+FiniteDiffWENO5.jl is a registered package and may be installed directly with the following command in the Julia REPL
+
+```julia-repl
+julia>]
+  pkg> add FiniteDiffWENO5
+  pkg> test FiniteDiffWENO5
+```

examples/1D/1D_linear.jl

@@ -1,13 +1,13 @@
 using FiniteDiffWENO5
 using GLMakie
 
-function main(nx=400)
+function 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)
+    x = range(x_min, stop = x_max, length = nx)
 
     # Courant number
     CFL = 0.4
@@ -16,20 +16,20 @@ function main(nx=400)
     # Parameters for Shu test
     z = -0.7
     δ = 0.005
-    β = log(2)/(36*δ^2)
+    β = 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))
+    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 + δ))
+    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)
@@ -42,18 +42,18 @@ function main(nx=400)
 
     # 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 + δ))
+    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)
+    weno = WENOScheme(u; boundary = (2, 2), stag = true, multithreading = false)
 
     # advection velocity
-    a = (;x=ones(nx+1))
+    a = (; x = ones(nx + 1))
 
     # grid size
     Δx = x[2] - x[1]
-    Δt = CFL*Δx^(5/3)
+    Δt = CFL * Δx^(5 / 3)
 
     tmax = period * (Lx + Δx) / maximum(abs.(a.x))
 
@@ -69,14 +69,14 @@ function main(nx=400)
         end
     end
 
-    f = Figure(size = (800, 600), dpi=400)
+    f = Figure(size = (800, 600), dpi = 400)
     ax = Axis(f[1, 1], title = "1D linear advection after $period periods", xlabel = "x", ylabel = "u")
     lines!(ax, x, u0_vec, label = "Exact", linestyle = :dash, color = :red)
     scatter!(ax, x, u, label = "WENO5")
     xlims!(ax, x_min, x_max)
     axislegend(ax)
     save("1D_linear_advection.png", f)
-    display(f)
+    return display(f)
 end
 
-main()
+main()

examples/1D/1D_linear_chmy.jl

@@ -3,17 +3,17 @@ using CairoMakie
 using Chmy
 using KernelAbstractions
 
-function main(backend=CPU(), nx=400)
+function main(backend = CPU(), nx = 400)
 
     arch = Arch(backend)
 
     x_min = -1.0
     x_max = 1.0
     Lx = x_max - x_min
 
-    x = range(x_min, stop=x_max, length=nx)
+    x = range(x_min, stop = x_max, length = nx)
 
-    grid   = UniformGrid(arch; origin=(x_min,), extent=(Lx,), dims=(nx,))
+    grid = UniformGrid(arch; origin = (x_min,), extent = (Lx,), dims = (nx,))
 
     # Courant number
     CFL = 0.7
@@ -22,20 +22,20 @@ function main(backend=CPU(), nx=400)
     # Parameters
     z = -0.7
     δ = 0.005
-    β = log(2)/(36*δ^2)
+    β = 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))
+    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 + δ))
+    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)
@@ -48,23 +48,23 @@ function main(backend=CPU(), nx=400)
 
     # 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 + δ))
+    u0_vec[idx] .= (1 / 6) .* (F(x[idx], α, a - δ) .+ 4 .* F(x[idx], α, a) .+ F(x[idx], α, a + δ))
 
 
     u = Field(backend, grid, Center())
     set!(u, u0_vec)
-    weno = WENOScheme(u, grid; boundary=(2, 2), stag=true)
+    weno = WENOScheme(u, grid; boundary = (2, 2), stag = true)
 
     # advection velocity
-    a_vec = ones(nx+1) .* -1
+    a_vec = ones(nx + 1) .* -1
     a = VectorField(backend, grid)
     set!(a, a_vec)
 
     # grid size
     Δx = x[2] - x[1]
-    Δt = CFL * Δx^(5/3)
+    Δt = CFL * Δx^(5 / 3)
 
-    tmax = period * (Lx+Δx) / maximum(abs.(a.x))
+    tmax = period * (Lx + Δx) / maximum(abs.(a.x))
 
     t = 0
 
@@ -86,7 +86,7 @@ function main(backend=CPU(), nx=400)
     scatter!(ax, x, interior(u) |> Array, label = "WENO5")
     xlims!(ax, x_min, x_max)
     axislegend(ax)
-    display(f)
+    return display(f)
 end
 
-main()
+main()