diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 059e10e..acf0628 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -32,11 +32,11 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v2 + - uses: julia-actions/cache@v3 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 diff --git a/.github/workflows/Runic.yml b/.github/workflows/Runic.yml index 806ee67..3ee92af 100644 --- a/.github/workflows/Runic.yml +++ b/.github/workflows/Runic.yml @@ -12,7 +12,7 @@ jobs: name: Runic runs-on: ubuntu-latest steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 - uses: fredrikekre/runic-action@v1 with: version: '1' \ No newline at end of file diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index fdcf679..53a5815 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -19,7 +19,7 @@ jobs: statuses: write runs-on: ubuntu-latest steps: - - uses: actions/checkout@v5 + - uses: actions/checkout@v6 - uses: julia-actions/setup-julia@latest with: version: '1' diff --git a/Project.toml b/Project.toml index b1e8176..9c2f200 100644 --- a/Project.toml +++ b/Project.toml @@ -1,23 +1,21 @@ name = "FiniteDiffWENO5" uuid = "4b1b4241-afed-47c3-a843-c6e03821b375" -version = "0.0.2" +version = "0.0.5" authors = ["Hugo Dominguez "] [deps] MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" -UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [weakdeps] -KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9" +KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" [extensions] -KAExt = ["KernelAbstractions"] ChmyExt = ["KernelAbstractions", "Chmy"] +KAExt = ["KernelAbstractions"] [compat] Chmy = "0.1" KernelAbstractions = "0.9" MuladdMacro = "0.2" -UnPack = "1.0" julia = "1.10" diff --git a/README.md b/README.md index d8a8475..877573c 100644 --- a/README.md +++ b/README.md @@ -14,6 +14,22 @@ Currently, the package focuses on non-conservative form of the advection terms ( 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 2 separate extensions. +## Installation + +FiniteDiffWENO5.jl is a registered Julia package and can be installed directly using the package manager: + +```julia-repl +julia>] + pkg> add FiniteDiffWENO5 +``` + +And you can test the package with: + +```julia-repl +julia>] + pkg> test FiniteDiffWENO5 +``` + ## Features The package currently exports only two main functions: `WENOScheme()`, that is used to create a WENO scheme struct containing all the necessary information for the WENO method, and `WENO_step!()`, that performs one step of the time integration using the WENO-Z method and a 3rd-order Runge-Kutta method. The grid and the initial condition must be defined by the user. @@ -108,3 +124,51 @@ end Which produces the following result: ![](/docs/src/assets/1D_linear_advection.png) + +## Multi-field advection + +If you have multiple scalar fields (e.g. different chemical components) that share the same velocity, you can advect them all in a single call by passing a tuple of arrays. The same `WENOScheme` buffers are reused for each field, so there is no extra memory overhead. Each field gets its own `u_min` / `u_max` bounds for the Zhang-Shu limiter. + +```julia +using FiniteDiffWENO5 + +nx, ny = 200, 200 +Lx = 1.0 +Δx, Δy = Lx / nx, Lx / ny + +# Three chemical components with different initial distributions +c1 = rand(nx, ny) +c2 = zeros(nx, ny) +c3 = ones(nx, ny) + +# Shared velocity field +v = (; x = ones(nx, ny), y = 0.5 .* ones(nx, ny)) + +# Create the WENO scheme from any one of the fields (they must all have the same size and type) +weno = WENOScheme(c1; boundary = (2, 2, 2, 2), stag = false) + +Δt = 0.7 * min(Δx, Δy)^(5 / 3) + +# Advect all three fields in one call with per-field limiter bounds +WENO_step!((c1, c2, c3), v, weno, Δt, Δx, Δy; + u_min = (0.0, 0.0, 0.0), + u_max = (1.0, 1.0, 1.0)) +``` + +This also works with the KernelAbstractions and Chmy backends: + +```julia +# KernelAbstractions +WENO_step!((c1, c2, c3), v, weno, Δt, Δx, Δy, backend; + u_min = (0.0, 0.0, 0.0), u_max = (1.0, 1.0, 1.0)) + +# Chmy +WENO_step!((c1, c2, c3), v, weno, Δt, Δx, Δy, grid, arch; + u_min = (0.0, 0.0, 0.0), u_max = (1.0, 1.0, 1.0)) +``` + +## Funding & author + +The development of this package was supported by the TRIGGER project funded by the German Federal Ministry for Economic Affairs and Energy (BMWK). + +Author: Hugo Dominguez (hdomingu@uni-mainz.de). diff --git a/docs/src/GPU.md b/docs/src/GPU.md new file mode 100644 index 0000000..5d33ac9 --- /dev/null +++ b/docs/src/GPU.md @@ -0,0 +1,106 @@ + +## Running on GPU + +To use this package on a GPU, you need to use [KernelAbstractions.jl](https://github.com/JuliaGPU/KernelAbstractions.jl) to define the arrays and run the computations. +The KernelAbstractions.jl package introduces a macro-based programming model that simplifies vendor-specific GPU programming by abstracting away its complexities. This allows hardware-independent kernels to be written that can be compiled and executed on different device backends without altering the high-level code or compromising performance. + + +This will load an extension from FiniteDiffWENO5GPU.jl: + +```julia +using KernelAbstractions +using FiniteDiffWENO5 +``` + +That way, the two functions `WENOScheme()` and `WENO_step!()` will automatically dispatch to the array types if the argument `backend` is provided. See the KernelAbstractions.jl documentation for more details on how to set up the GPU backend. + + + + +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) \ No newline at end of file diff --git a/docs/src/GettingStarted.md b/docs/src/GettingStarted.md index 03ca2fe..b49c1a6 100644 --- a/docs/src/GettingStarted.md +++ b/docs/src/GettingStarted.md @@ -96,4 +96,50 @@ display(f) which outputs: -![1D advection](./assets/1D_linear_advection.png) \ No newline at end of file +![1D advection](./assets/1D_linear_advection.png) + +## Advecting multiple fields + +If you have multiple scalar fields (e.g. different chemical components) that share the same velocity, you can advect them all in a single call by passing a tuple of arrays. The same `WENOScheme` buffers are reused for each field, so there is no extra memory overhead. + +Each field gets its own `u_min` / `u_max` bounds for the Zhang-Shu limiter, passed as tuples matching the number of fields. + +```julia +using FiniteDiffWENO5 + +nx = 200 +ny = 200 +Lx = 1.0 +Δx = Lx / nx +Δy = Lx / ny + +# Three chemical components with different initial distributions +c1 = rand(nx, ny) # component 1 +c2 = zeros(nx, ny) # component 2 +c3 = ones(nx, ny) # component 3 + +# Shared velocity field +v = (; x = ones(nx, ny), y = 0.5 .* ones(nx, ny)) + +# Create the WENO scheme from any one of the fields (they must all have the same size/type) +weno = WENOScheme(c1; boundary = (2, 2, 2, 2), stag = false) + +Δt = 0.7 * min(Δx, Δy)^(5 / 3) + +# Advect all three fields in one call +WENO_step!((c1, c2, c3), v, weno, Δt, Δx, Δy; + u_min = (0.0, 0.0, 0.0), + u_max = (1.0, 1.0, 1.0)) +``` + +This also works with KernelAbstractions and Chmy backends — just pass the extra backend/grid arguments as usual: + +```julia +# KernelAbstractions +WENO_step!((c1, c2, c3), v, weno, Δt, Δx, Δy, backend; + u_min = (0.0, 0.0, 0.0), u_max = (1.0, 1.0, 1.0)) + +# Chmy +WENO_step!((c1, c2, c3), v, weno, Δt, Δx, Δy, grid, arch; + u_min = (0.0, 0.0, 0.0), u_max = (1.0, 1.0, 1.0)) +``` \ No newline at end of file diff --git a/examples/2D/2D_multi_field.jl b/examples/2D/2D_multi_field.jl new file mode 100644 index 0000000..b1b948c --- /dev/null +++ b/examples/2D/2D_multi_field.jl @@ -0,0 +1,92 @@ +using FiniteDiffWENO5 +using GLMakie + +function main(; nx = 400, ny = 400) + + Lx = 1.0 + Δx = Lx / nx + Δy = Lx / ny + + # Courant number + CFL = 0.7 + period = 1 + + # Grid + x = range(0, length = nx, stop = Lx) + y = range(0, length = ny, stop = Lx) + grid = (x .* ones(ny)', ones(nx) .* y') + + # Shared velocity field + vx0 = ones(nx, ny) + vy0 = ones(nx, ny) + v = (; x = vy0, y = vx0) + + # Three chemical components with different initial conditions + x0 = 1 / 4 + c_width1 = 0.08 + c_width2 = 0.06 + c_width3 = 0.1 + + c1 = zeros(ny, nx) + c2 = zeros(ny, nx) + c3 = zeros(ny, nx) + + for I in CartesianIndices((ny, nx)) + c1[I] = sign(exp(-((grid[1][I] - x0)^2 + (grid[2][I]' - x0)^2) / c_width1^2) - 0.5) * 0.5 + 0.5 + c2[I] = exp(-((grid[1][I] - 0.5)^2 + (grid[2][I]' - 0.5)^2) / c_width2^2) + c3[I] = exp(-((grid[1][I] - 0.75)^2 + (grid[2][I]' - 0.25)^2) / c_width3^2) + end + + c1_0 = copy(c1) + c2_0 = copy(c2) + c3_0 = copy(c3) + + # Create a single WENOScheme shared by all fields + weno = WENOScheme(c1; boundary = (2, 2, 2, 2), stag = false, multithreading = true) + + Δt = CFL * min(Δx, Δy)^(5 / 3) + tmax = period * Lx / max(maximum(abs.(vx0)), maximum(abs.(vy0))) + + t = 0 + counter = 0 + + # Plot setup + f = Figure(size = (1200, 400)) + ax1 = Axis(f[1, 1], title = "Component 1") + ax2 = Axis(f[1, 2], title = "Component 2") + ax3 = Axis(f[1, 3], title = "Component 3") + c1_obs = Observable(c1_0) + c2_obs = Observable(c2_0) + c3_obs = Observable(c3_0) + heatmap!(ax1, x, y, c1_obs; colormap = cgrad(:roma, rev = true), colorrange = (0, 1)) + heatmap!(ax2, x, y, c2_obs; colormap = cgrad(:roma, rev = true), colorrange = (0, 1)) + heatmap!(ax3, x, y, c3_obs; colormap = cgrad(:roma, rev = true), colorrange = (0, 1)) + display(f) + + while t < tmax + # Advect all 3 components in a single call + WENO_step!( + (c1, c2, c3), v, weno, Δt, Δx, Δy; + u_min = (0.0, 0.0, 0.0), + u_max = (1.0, 1.0, 1.0) + ) + + t += Δt + + if t + Δt > tmax + Δt = tmax - t + end + + if counter % 100 == 0 + c1_obs[] = c1 + c2_obs[] = c2 + c3_obs[] = c3 + end + + counter += 1 + end + + return nothing +end + +main(nx = 400, ny = 400) diff --git a/ext/ChmyExt.jl b/ext/ChmyExt.jl index 3373108..dbf7c20 100644 --- a/ext/ChmyExt.jl +++ b/ext/ChmyExt.jl @@ -2,7 +2,6 @@ module ChmyExt using FiniteDiffWENO5 using FiniteDiffWENO5: zhang_shu_limit, weno5_reconstruction_upwind, weno5_reconstruction_downwind using MuladdMacro -using UnPack using Chmy using KernelAbstractions @@ -103,7 +102,7 @@ function WENO_step!(u::T_field, v::NamedTuple{(:x,), <:Tuple{<:AbstractField{<:R nx = grid.axes[1].length Δx_ = inv(Δx) - @unpack ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode = weno + (; ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode) = weno if !upwind_mode @@ -164,7 +163,7 @@ function WENO_step!(u::T_field, v::NamedTuple{(:x, :y), <:Tuple{Vararg{AbstractF Δx_ = inv(Δx) Δy_ = inv(Δy) - @unpack ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode = weno + (; ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode) = weno if !upwind_mode launch(arch, grid, WENO_flux_KA_2D_x => (fl.x, fr.x, u, boundary, nx, χ, γ, ζ, ϵ, lim_ZS, u_min, u_max, grid)) @@ -232,7 +231,7 @@ function WENO_step!(u::T_field, v::NamedTuple{(:x, :y, :z), <:Tuple{Vararg{Abstr Δy_ = inv(Δy) Δz_ = inv(Δz) - @unpack ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode = weno + (; ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode) = weno if !upwind_mode launch(arch, grid, WENO_flux_KA_3D_x => (fl.x, fr.x, u, boundary, nx, χ, γ, ζ, ϵ, lim_ZS, u_min, u_max, grid)) @@ -263,4 +262,43 @@ function WENO_step!(u::T_field, v::NamedTuple{(:x, :y, :z), <:Tuple{Vararg{Abstr end +# ── Multi-field overloads ────────────────────────────────────────────── + +""" + WENO_step!(u::NTuple{NF}, v, weno, Δt, Δx, grid, arch; u_min, u_max) + +Advance multiple 1D Chmy fields sharing the same velocity and WENOScheme buffers. +""" +function WENO_step!(u::Tuple{Vararg{<:AbstractField{<:Real, 1}}}, v::NamedTuple{(:x,), <:Tuple{<:AbstractField{<:Real, 1}}}, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, grid::StructuredGrid, arch; u_min::Tuple{Vararg{Real}}, u_max::Tuple{Vararg{Real}}) + for i in eachindex(u) + WENO_step!(u[i], v, weno, Δt, Δx, grid, arch; u_min = u_min[i], u_max = u_max[i]) + end + return nothing +end + +""" + WENO_step!(u::NTuple{NF}, v, weno, Δt, Δx, Δy, grid, arch; u_min, u_max) + +Advance multiple 2D Chmy fields sharing the same velocity and WENOScheme buffers. +""" +function WENO_step!(u::Tuple{Vararg{<:AbstractField{<:Real, 2}}}, v::NamedTuple{(:x, :y), <:Tuple{Vararg{AbstractField{<:Real}, 2}}}, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, Δy, grid::StructuredGrid, arch; u_min::Tuple{Vararg{Real}}, u_max::Tuple{Vararg{Real}}) + for i in eachindex(u) + WENO_step!(u[i], v, weno, Δt, Δx, Δy, grid, arch; u_min = u_min[i], u_max = u_max[i]) + end + return nothing +end + +""" + WENO_step!(u::NTuple{NF}, v, weno, Δt, Δx, Δy, Δz, grid, arch; u_min, u_max) + +Advance multiple 3D Chmy fields sharing the same velocity and WENOScheme buffers. +""" +function WENO_step!(u::Tuple{Vararg{<:AbstractArray{<:Real, 3}}}, v::NamedTuple{(:x, :y, :z), <:Tuple{Vararg{AbstractField{<:Real}, 3}}}, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, Δy, Δz, grid::StructuredGrid, arch; u_min::Tuple{Vararg{Real}}, u_max::Tuple{Vararg{Real}}) + for i in eachindex(u) + WENO_step!(u[i], v, weno, Δt, Δx, Δy, Δz, grid, arch; u_min = u_min[i], u_max = u_max[i]) + end + return nothing +end + + end diff --git a/ext/KAExt.jl b/ext/KAExt.jl index 4f33b85..8751551 100644 --- a/ext/KAExt.jl +++ b/ext/KAExt.jl @@ -2,7 +2,6 @@ module KAExt using FiniteDiffWENO5 using FiniteDiffWENO5: zhang_shu_limit, weno5_reconstruction_upwind, weno5_reconstruction_downwind using MuladdMacro -using UnPack using KernelAbstractions import FiniteDiffWENO5: WENOScheme, WENO_step! @@ -99,7 +98,7 @@ function WENO_step!(u::T_KA, v::NamedTuple{(:x,), <:Tuple{<:AbstractArray{<:Real #! do things here for halos and such for clusters for boundaries probably - @unpack ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode = weno + (; ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode) = weno nx = size(u, 1) Δx_ = inv(Δx) @@ -163,7 +162,7 @@ function WENO_step!(u::T_KA, v::NamedTuple{(:x, :y), <:Tuple{Vararg{AbstractArra #! do things here for halos and such for clusters for boundaries probably - @unpack ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode = weno + (; ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode) = weno nx = size(u, 1) ny = size(u, 2) @@ -245,7 +244,7 @@ function WENO_step!(u::T_KA, v::NamedTuple{(:x, :y, :z), <:Tuple{Vararg{Abstract Δy_ = inv(Δy) Δz_ = inv(Δz) - @unpack ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode = weno + (; ut, du, fl, fr, stag, lim_ZS, boundary, χ, γ, ζ, ϵ, upwind_mode) = weno if !upwind_mode @@ -293,4 +292,43 @@ function WENO_step!(u::T_KA, v::NamedTuple{(:x, :y, :z), <:Tuple{Vararg{Abstract end +# ── Multi-field overloads ────────────────────────────────────────────── + +""" + WENO_step!(u::NTuple{NF}, v, weno, Δt, Δx, backend; u_min, u_max) + +Advance multiple 1D fields sharing the same velocity and WENOScheme buffers. +""" +function WENO_step!(u::Tuple{Vararg{<:AbstractVector{<:Real}}}, v::NamedTuple{(:x,), <:Tuple{<:AbstractArray{<:Real}}}, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, backend::Backend; u_min::Tuple{Vararg{Real}}, u_max::Tuple{Vararg{Real}}) + for i in eachindex(u) + WENO_step!(u[i], v, weno, Δt, Δx, backend; u_min = u_min[i], u_max = u_max[i]) + end + return nothing +end + +""" + WENO_step!(u::NTuple{NF}, v, weno, Δt, Δx, Δy, backend; u_min, u_max) + +Advance multiple 2D fields sharing the same velocity and WENOScheme buffers. +""" +function WENO_step!(u::Tuple{Vararg{<:AbstractArray{<:Real, 2}}}, v::NamedTuple{(:x, :y), <:Tuple{Vararg{AbstractArray{<:Real}, 2}}}, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, Δy, backend::Backend; u_min::Tuple{Vararg{Real}}, u_max::Tuple{Vararg{Real}}) + for i in eachindex(u) + WENO_step!(u[i], v, weno, Δt, Δx, Δy, backend; u_min = u_min[i], u_max = u_max[i]) + end + return nothing +end + +""" + WENO_step!(u::NTuple{NF}, v, weno, Δt, Δx, Δy, Δz, backend; u_min, u_max) + +Advance multiple 3D fields sharing the same velocity and WENOScheme buffers. +""" +function WENO_step!(u::Tuple{Vararg{<:AbstractArray{<:Real, 3}}}, v::NamedTuple{(:x, :y, :z), <:Tuple{Vararg{AbstractArray{<:Real}, 3}}}, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, Δy, Δz, backend::Backend; u_min::Tuple{Vararg{Real}}, u_max::Tuple{Vararg{Real}}) + for i in eachindex(u) + WENO_step!(u[i], v, weno, Δt, Δx, Δy, Δz, backend; u_min = u_min[i], u_max = u_max[i]) + end + return nothing +end + + end diff --git a/src/1D/semi_discretisation_1D.jl b/src/1D/semi_discretisation_1D.jl index 7cf6a13..cdeac10 100644 --- a/src/1D/semi_discretisation_1D.jl +++ b/src/1D/semi_discretisation_1D.jl @@ -1,5 +1,5 @@ function WENO_flux!(fl, fr, u, weno, nx, u_min, u_max) - @unpack boundary, χ, γ, ζ, ϵ, multithreading, lim_ZS = weno + (; boundary, χ, γ, ζ, ϵ, multithreading, lim_ZS) = weno bL = Val(boundary[1]) bR = Val(boundary[2]) @@ -37,7 +37,7 @@ end function semi_discretisation_weno5!(du::T, v, weno::WENOScheme, Δx_) where {T <: AbstractArray{<:Real, 1}} - @unpack fl, fr, stag, multithreading = weno + (; fl, fr, stag, multithreading) = weno # use staggered grid or not for the velocities if stag @@ -70,7 +70,7 @@ using velocity field `v` on a staggered or collocated grid. - Uses the boundary conditions stored in `weno.boundary`. """ function upwind_update_1D!(u, v, weno, nx, Δx_, Δt) - @unpack boundary, stag, multithreading = weno + (; boundary, stag, multithreading) = weno bL = Val(boundary[1]) bR = Val(boundary[2]) diff --git a/src/1D/time_stepping.jl b/src/1D/time_stepping.jl index 17a47e4..3c44e51 100644 --- a/src/1D/time_stepping.jl +++ b/src/1D/time_stepping.jl @@ -24,7 +24,7 @@ function WENO_step!(u::T, v::NamedTuple{(:x,), <:Tuple{<:Vector{<:Real}}}, weno: nx = size(u, 1) Δx_ = inv(Δx) - @unpack ut, du, stag, fl, fr, multithreading, upwind_mode = weno + (; ut, du, stag, fl, fr, multithreading, upwind_mode) = weno if !upwind_mode WENO_flux!(fl, fr, u, weno, nx, u_min, u_max) @@ -55,3 +55,21 @@ function WENO_step!(u::T, v::NamedTuple{(:x,), <:Tuple{<:Vector{<:Real}}}, weno: return nothing end + +""" + WENO_step!(u::Tuple{Vararg{Vector{<:Real}}}, + v::NamedTuple{(:x,), <:Tuple{<:Vector{<:Real}}}, + weno::WENOScheme, + Δt, Δx; + u_min::Tuple{Vararg{Real}}, + u_max::Tuple{Vararg{Real}}) + +Advance multiple fields `u = (c1, c2, ...)` by one time step, all sharing the same velocity `v` and `WENOScheme` buffers. +Each field is advected sequentially with its own `u_min` / `u_max` bounds for the Zhang-Shu limiter. +""" +function WENO_step!(u::Tuple{Vararg{Vector{<:Real}}}, v::NamedTuple{(:x,), <:Tuple{<:Vector{<:Real}}}, weno::WENOScheme, Δt, Δx; u_min::Tuple{Vararg{Real}}, u_max::Tuple{Vararg{Real}}) + for i in eachindex(u) + WENO_step!(u[i], v, weno, Δt, Δx; u_min = u_min[i], u_max = u_max[i]) + end + return nothing +end diff --git a/src/2D/semi_discretisation_2D.jl b/src/2D/semi_discretisation_2D.jl index 7807895..3a34422 100644 --- a/src/2D/semi_discretisation_2D.jl +++ b/src/2D/semi_discretisation_2D.jl @@ -1,5 +1,5 @@ function WENO_flux!(fl, fr, u, weno, nx, ny, u_max, u_min) - @unpack boundary, χ, γ, ζ, ϵ, multithreading, lim_ZS = weno + (; boundary, χ, γ, ζ, ϵ, multithreading, lim_ZS) = weno bLx = Val(boundary[1]) bRx = Val(boundary[2]) @@ -93,7 +93,7 @@ end function semi_discretisation_weno5!(du::T, v, weno::WENOScheme, Δx_, Δy_) where {T <: AbstractArray{<:Real, 2}} - @unpack fl, fr, stag, multithreading = weno + (; fl, fr, stag, multithreading) = weno # use staggered grid or not for the velocities if stag @@ -132,7 +132,7 @@ end function upwind_update_2D!( u, v, weno, nx, ny, Δx_, Δy_, Δt ) - @unpack boundary, stag, multithreading = weno + (; boundary, stag, multithreading) = weno bLx = Val(boundary[1]) bRx = Val(boundary[2]) diff --git a/src/2D/time_stepping.jl b/src/2D/time_stepping.jl index 607530d..b604bcb 100644 --- a/src/2D/time_stepping.jl +++ b/src/2D/time_stepping.jl @@ -25,7 +25,7 @@ function WENO_step!(u::T, v::NamedTuple{(:x, :y), <:Tuple{Vararg{AbstractArray{< nx, ny = size(u, 1), size(u, 2) Δx_, Δy_ = inv(Δx), inv(Δy) - @unpack ut, du, stag, fl, fr, multithreading, upwind_mode = weno + (; ut, du, stag, fl, fr, multithreading, upwind_mode) = weno if !upwind_mode WENO_flux!(fl, fr, u, weno, nx, ny, u_max, u_min) @@ -55,3 +55,21 @@ function WENO_step!(u::T, v::NamedTuple{(:x, :y), <:Tuple{Vararg{AbstractArray{< return nothing end + +""" + WENO_step!(u::Tuple{Vararg{Array{<:Real, 2}}}, + v::NamedTuple{(:x, :y), <:Tuple{Vararg{AbstractArray{<:Real}, 2}}}, + weno::WENOScheme, + Δt, Δx, Δy; + u_min::Tuple{Vararg{Real}}, + u_max::Tuple{Vararg{Real}}) + +Advance multiple fields `u = (c1, c2, ...)` by one time step, all sharing the same velocity `v` and `WENOScheme` buffers. +Each field is advected sequentially with its own `u_min` / `u_max` bounds for the Zhang-Shu limiter. +""" +function WENO_step!(u::Tuple{Vararg{Array{<:Real, 2}}}, v::NamedTuple{(:x, :y), <:Tuple{Vararg{AbstractArray{<:Real}, 2}}}, weno::WENOScheme, Δt, Δx, Δy; u_min::Tuple{Vararg{Real}}, u_max::Tuple{Vararg{Real}}) + for i in eachindex(u) + WENO_step!(u[i], v, weno, Δt, Δx, Δy; u_min = u_min[i], u_max = u_max[i]) + end + return nothing +end diff --git a/src/3D/semi_discretisation_3D.jl b/src/3D/semi_discretisation_3D.jl index 97632e4..ede9a38 100644 --- a/src/3D/semi_discretisation_3D.jl +++ b/src/3D/semi_discretisation_3D.jl @@ -1,5 +1,5 @@ function WENO_flux!(fl, fr, u, weno, nx, ny, nz, u_min, u_max) - @unpack boundary, χ, γ, ζ, ϵ, multithreading, lim_ZS = weno + (; boundary, χ, γ, ζ, ϵ, multithreading, lim_ZS) = weno bLx = Val(boundary[1]) bRx = Val(boundary[2]) @@ -149,7 +149,7 @@ end function semi_discretisation_weno5!(du::T, v, weno::WENOScheme, Δx_, Δy_, Δz_) where {T <: AbstractArray{<:Real, 3}} - @unpack fl, fr, stag, multithreading = weno + (; fl, fr, stag, multithreading) = weno # use staggered grid or not for the velocities if stag @@ -197,7 +197,7 @@ end function upwind_update_3D!( u, v, weno, nx, ny, nz, Δx_, Δy_, Δz_, Δt ) - @unpack boundary, stag, multithreading = weno + (; boundary, stag, multithreading) = weno bLx = Val(boundary[1]) bRx = Val(boundary[2]) diff --git a/src/3D/time_stepping.jl b/src/3D/time_stepping.jl index ad90e15..938f2d0 100644 --- a/src/3D/time_stepping.jl +++ b/src/3D/time_stepping.jl @@ -26,7 +26,7 @@ function WENO_step!(u::T, v::NamedTuple{(:x, :y, :z), <:Tuple{Vararg{Array{<:Rea nx, ny, nz = size(u, 1), size(u, 2), size(u, 3) Δx_, Δy_, Δz_ = inv(Δx), inv(Δy), inv(Δz) - @unpack ut, du, stag, fl, fr, multithreading, upwind_mode = weno + (; ut, du, stag, fl, fr, multithreading, upwind_mode) = weno if !upwind_mode WENO_flux!(fl, fr, u, weno, nx, ny, nz, u_min, u_max) @@ -56,3 +56,21 @@ function WENO_step!(u::T, v::NamedTuple{(:x, :y, :z), <:Tuple{Vararg{Array{<:Rea return nothing end + +""" + WENO_step!(u::Tuple{Vararg{Array{<:Real, 3}}}, + v::NamedTuple{(:x, :y, :z), <:Tuple{Vararg{Array{<:Real}, 3}}}, + weno::WENOScheme, + Δt, Δx, Δy, Δz; + u_min::Tuple{Vararg{Real}}, + u_max::Tuple{Vararg{Real}}) + +Advance multiple fields `u = (c1, c2, ...)` by one time step, all sharing the same velocity `v` and `WENOScheme` buffers. +Each field is advected sequentially with its own `u_min` / `u_max` bounds for the Zhang-Shu limiter. +""" +function WENO_step!(u::Tuple{Vararg{Array{<:Real, 3}}}, v::NamedTuple{(:x, :y, :z), <:Tuple{Vararg{Array{<:Real}, 3}}}, weno::WENOScheme, Δt, Δx, Δy, Δz; u_min::Tuple{Vararg{Real}}, u_max::Tuple{Vararg{Real}}) + for i in eachindex(u) + WENO_step!(u[i], v, weno, Δt, Δx, Δy, Δz; u_min = u_min[i], u_max = u_max[i]) + end + return nothing +end diff --git a/src/FiniteDiffWENO5.jl b/src/FiniteDiffWENO5.jl index 55e3e55..4c33366 100644 --- a/src/FiniteDiffWENO5.jl +++ b/src/FiniteDiffWENO5.jl @@ -1,6 +1,5 @@ module FiniteDiffWENO5 -using UnPack using MuladdMacro export WENOScheme, WENO_step! diff --git a/test/Project.toml b/test/Project.toml index e6d34b6..9536f83 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,7 +4,6 @@ KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" [compat] Aqua = "0.8" \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 4a7b560..c94b9db 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,7 +2,6 @@ using Test using FiniteDiffWENO5 using KernelAbstractions using Chmy -using UnPack function runtests() files = readdir(@__DIR__) diff --git a/test/test_multi_field.jl b/test/test_multi_field.jl new file mode 100644 index 0000000..055cf5e --- /dev/null +++ b/test/test_multi_field.jl @@ -0,0 +1,296 @@ +@testset "Multi-field advection tests" begin + + @testset "1D multi-field equivalence" begin + nx = 200 + x_min = -1.0 + x_max = 1.0 + Lx = x_max - x_min + x = range(x_min, stop = x_max, length = nx) + + CFL = 0.4 + period = 1 + + # Two different initial conditions + u1_init = @. exp(-((x - 0.0)^2) / 0.01) + u2_init = @. exp(-((x + 0.3)^2) / 0.02) + u3_init = @. exp(-((x - 0.4)^2) / 0.005) + + # Single-field references + u1_single = copy(u1_init) + u2_single = copy(u2_init) + u3_single = copy(u3_init) + + # Multi-field tuple + u1_multi = copy(u1_init) + u2_multi = copy(u2_init) + u3_multi = copy(u3_init) + + weno = WENOScheme(u1_init; boundary = (2, 2), stag = true) + + a = (; x = ones(nx + 1)) + Δx = x[2] - x[1] + Δt = CFL * Δx^(5 / 3) + tmax = period * (Lx + Δx) / maximum(abs.(a.x)) + + t = 0.0 + while t < tmax + # Single-field calls + WENO_step!(u1_single, a, weno, Δt, Δx; u_min = 0.0, u_max = 1.0) + WENO_step!(u2_single, a, weno, Δt, Δx; u_min = 0.0, u_max = 1.0) + WENO_step!(u3_single, a, weno, Δt, Δx; u_min = -1.0, u_max = 2.0) + + # Multi-field call + WENO_step!( + (u1_multi, u2_multi, u3_multi), a, weno, Δt, Δx; + u_min = (0.0, 0.0, -1.0), u_max = (1.0, 1.0, 2.0) + ) + + t += Δt + if t + Δt > tmax + Δt = tmax - t + end + end + + @test u1_multi ≈ u1_single + @test u2_multi ≈ u2_single + @test u3_multi ≈ u3_single + end + + + @testset "2D multi-field equivalence" begin + nx = 50 + ny = 50 + Lx = 1.0 + Δx = Lx / nx + Δy = Lx / ny + + x = range(0, length = nx, stop = Lx) + y = range(0, length = ny, stop = Lx) + grid_array = (x .* ones(ny)', ones(nx) .* y') + + CFL = 0.7 + period = 1 + + x0 = 1 / 4 + c1 = 0.08 + c2 = 0.12 + + u1_init = zeros(ny, nx) + u2_init = zeros(ny, nx) + for I in CartesianIndices((ny, nx)) + u1_init[I] = exp(-((grid_array[1][I] - x0)^2 + (grid_array[2][I]' - x0)^2) / c1^2) + u2_init[I] = exp(-((grid_array[1][I] - 0.5)^2 + (grid_array[2][I]' - 0.5)^2) / c2^2) + end + + vx0 = ones(nx, ny) + vy0 = ones(nx, ny) + v = (; x = vx0, y = vy0) + + # Single-field references + u1_single = copy(u1_init) + u2_single = copy(u2_init) + + # Multi-field + u1_multi = copy(u1_init) + u2_multi = copy(u2_init) + + weno = WENOScheme(u1_init; boundary = (2, 2, 2, 2), stag = false, multithreading = false) + + Δt = CFL * min(Δx, Δy)^(5 / 3) + tmax = period * Lx / max(maximum(abs.(vx0)), maximum(abs.(vy0))) + + t = 0.0 + while t < tmax + WENO_step!(u1_single, v, weno, Δt, Δx, Δy; u_min = 0.0, u_max = 1.0) + WENO_step!(u2_single, v, weno, Δt, Δx, Δy; u_min = 0.0, u_max = 1.0) + + WENO_step!( + (u1_multi, u2_multi), v, weno, Δt, Δx, Δy; + u_min = (0.0, 0.0), u_max = (1.0, 1.0) + ) + + t += Δt + if t + Δt > tmax + Δt = tmax - t + end + end + + @test u1_multi ≈ u1_single + @test u2_multi ≈ u2_single + end + + + @testset "3D multi-field equivalence" begin + nx = 20 + ny = 20 + nz = 20 + + L = 1.0 + Δx = L / nx + Δy = L / ny + Δz = L / nz + + x = range(0, stop = L, length = nx) + y = range(0, stop = L, length = ny) + z = range(0, stop = L, length = nz) + + X = reshape(x, 1, 1, nx) .* ones(ny, nz, 1) + Y = reshape(y, 1, ny, 1) .* ones(nx, 1, nz) + Z = reshape(z, nz, 1, 1) .* ones(1, ny, nx) + + CFL = 0.7 + period = 1 + + u1_init = zeros(ny, nx, nz) + u2_init = zeros(ny, nx, nz) + for I in CartesianIndices((ny, nx, nz)) + u1_init[I] = exp(-((X[I] - 0.25)^2 + (Y[I] - 0.25)^2 + (Z[I] - 0.5)^2) / 0.01) + u2_init[I] = exp(-((X[I] - 0.5)^2 + (Y[I] - 0.5)^2 + (Z[I] - 0.5)^2) / 0.02) + end + + vx0 = ones(size(X)) + vy0 = ones(size(Y)) + vz0 = zeros(size(Z)) + v = (; x = vx0, y = vy0, z = vz0) + + u1_single = copy(u1_init) + u2_single = copy(u2_init) + u1_multi = copy(u1_init) + u2_multi = copy(u2_init) + + weno = WENOScheme(u1_init; boundary = (2, 2, 2, 2, 2, 2), stag = false, multithreading = false) + + Δt = CFL * min(Δx, Δy, Δz)^(5 / 3) + tmax = period * L / max(maximum(abs.(vx0)), maximum(abs.(vy0))) + + t = 0.0 + while t < tmax + WENO_step!(u1_single, v, weno, Δt, Δx, Δy, Δz; u_min = 0.0, u_max = 1.0) + WENO_step!(u2_single, v, weno, Δt, Δx, Δy, Δz; u_min = 0.0, u_max = 1.0) + + WENO_step!( + (u1_multi, u2_multi), v, weno, Δt, Δx, Δy, Δz; + u_min = (0.0, 0.0), u_max = (1.0, 1.0) + ) + + t += Δt + if t + Δt > tmax + Δt = tmax - t + end + end + + @test u1_multi ≈ u1_single + @test u2_multi ≈ u2_single + end + + + @testset "1D multi-field KA CPU equivalence" begin + backend = CPU() + nx = 200 + x_min = -1.0 + x_max = 1.0 + Lx = x_max - x_min + x = range(x_min, stop = x_max, length = nx) + + CFL = 0.4 + period = 1 + + u1_init = KernelAbstractions.zeros(backend, Float64, nx) + u2_init = KernelAbstractions.zeros(backend, Float64, nx) + for i in 1:nx + u1_init[i] = exp(-((x[i] - 0.0)^2) / 0.01) + u2_init[i] = exp(-((x[i] + 0.3)^2) / 0.02) + end + + u1_single = copy(u1_init) + u2_single = copy(u2_init) + u1_multi = copy(u1_init) + u2_multi = copy(u2_init) + + weno = WENOScheme(u1_init, backend; boundary = (2, 2), stag = true) + + a = (; x = KernelAbstractions.ones(backend, Float64, nx + 1)) + Δx = x[2] - x[1] + Δt = CFL * Δx^(5 / 3) + tmax = period * (Lx + Δx) / 1.0 + + t = 0.0 + while t < tmax + WENO_step!(u1_single, a, weno, Δt, Δx, backend; u_min = 0.0, u_max = 1.0) + WENO_step!(u2_single, a, weno, Δt, Δx, backend; u_min = 0.0, u_max = 1.0) + + WENO_step!( + (u1_multi, u2_multi), a, weno, Δt, Δx, backend; + u_min = (0.0, 0.0), u_max = (1.0, 1.0) + ) + + t += Δt + if t + Δt > tmax + Δt = tmax - t + end + end + + @test u1_multi ≈ u1_single + @test u2_multi ≈ u2_single + end + + @testset "1D multi-field Chmy CPU equivalence" begin + backend = CPU() + nx = 200 + x_min = -1.0 + x_max = 1.0 + Lx = x_max - x_min + x = range(x_min, stop = x_max, length = nx) + + arch = Arch(backend) + grid = UniformGrid(arch; origin = (x_min,), extent = (Lx,), dims = (nx,)) + + CFL = 0.4 + period = 1 + + u1_init = zeros(nx) + u2_init = zeros(nx) + for i in 1:nx + u1_init[i] = exp(-((x[i] - 0.0)^2) / 0.01) + u2_init[i] = exp(-((x[i] + 0.3)^2) / 0.02) + end + + u1_single = Field(backend, grid, Center()) + u2_single = Field(backend, grid, Center()) + u1_multi = Field(backend, grid, Center()) + u2_multi = Field(backend, grid, Center()) + set!(u1_single, u1_init) + set!(u2_single, u2_init) + set!(u1_multi, u1_init) + set!(u2_multi, u2_init) + + weno = WENOScheme(u1_single, grid; boundary = (2, 2), stag = true) + + a_vec = ones(nx + 1) + a = VectorField(backend, grid) + set!(a, a_vec) + + Δx = x[2] - x[1] + Δt = CFL * Δx^(5 / 3) + tmax = period * (Lx + Δx) / maximum(abs.(a.x)) + + t = 0.0 + while t < tmax + WENO_step!(u1_single, a, weno, Δt, Δx, grid, arch; u_min = 0.0, u_max = 1.0) + WENO_step!(u2_single, a, weno, Δt, Δx, grid, arch; u_min = 0.0, u_max = 1.0) + + WENO_step!( + (u1_multi, u2_multi), a, weno, Δt, Δx, grid, arch; + u_min = (0.0, 0.0), u_max = (1.0, 1.0) + ) + + t += Δt + if t + Δt > tmax + Δt = tmax - t + end + end + + @test interior(u1_multi) ≈ interior(u1_single) + @test interior(u2_multi) ≈ interior(u2_single) + end +end diff --git a/test/test_reconstruction.jl b/test/test_reconstruction.jl index 317cb3b..5d6b44e 100644 --- a/test/test_reconstruction.jl +++ b/test/test_reconstruction.jl @@ -4,7 +4,7 @@ @testset "Constant field" begin u = fill(3.14, 5) weno = WENOScheme(u) - @unpack χ, γ, ζ, ϵ = weno + (; χ, γ, ζ, ϵ) = weno f_up = FiniteDiffWENO5.weno5_reconstruction_upwind(u..., χ, γ, ζ, ϵ) f_down = FiniteDiffWENO5.weno5_reconstruction_downwind(u..., χ, γ, ζ, ϵ) @@ -17,7 +17,7 @@ @testset "Linear field" begin u = [1.0, 2.0, 3.0, 4.0, 5.0] weno = WENOScheme(u) - @unpack χ, γ, ζ, ϵ = weno + (; χ, γ, ζ, ϵ) = weno f_up = FiniteDiffWENO5.weno5_reconstruction_upwind(u..., χ, γ, ζ, ϵ) f_down = FiniteDiffWENO5.weno5_reconstruction_downwind(u..., χ, γ, ζ, ϵ) @@ -31,7 +31,7 @@ x = -2:2 u = float.(x .^ 2) # smooth convex profile weno = WENOScheme(float(u)) - @unpack χ, γ, ζ, ϵ = weno + (; χ, γ, ζ, ϵ) = weno fl = (2u[1] - 13u[2] + 47u[3] + 27u[4] - 3u[5]) / 60.0 fr = (-3u[1] + 27u[2] + 47u[3] - 13u[4] + 2u[5]) / 60.0 @@ -47,7 +47,7 @@ @testset "Discontinuous field" begin u = [1.0, 1.0, 1.0, 10.0, 10.0] weno = WENOScheme(u) - @unpack χ, γ, ζ, ϵ = weno + (; χ, γ, ζ, ϵ) = weno f_up = FiniteDiffWENO5.weno5_reconstruction_upwind(u..., χ, γ, ζ, ϵ) f_down = FiniteDiffWENO5.weno5_reconstruction_downwind(u..., χ, γ, ζ, ϵ)