add 3D

Author: Iddingsite

Project.toml

@@ -12,6 +12,7 @@ Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9"
 [extensions]
 ChmyExt1D = ["Chmy", "KernelAbstractions"]
 ChmyExt2D = ["Chmy", "KernelAbstractions"]
+ChmyExt3D = ["Chmy", "KernelAbstractions"]
 
 [compat]
 Chmy = "0.1"

examples/2D/2D_linear_chmy.jl

@@ -11,8 +11,6 @@ function main(;backend=CPU(), nx=400, ny=400, stag=true)
     Δx = Lx / nx
     Δy = Lx / ny
 
-    x = range(0, stop=Lx, length=nx)
-
     grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(Lx, Lx), dims=(nx, ny))
 
     # Courant number

examples/2D/2D_linear_rotation_chmy.jl

@@ -4,18 +4,13 @@ using KernelAbstractions
 using Plots
 
 function main(;backend=CPU(), nx=400, ny=400)
-# backend=CPU()
-# nx=400
-# ny=400
 
     arch = Arch(backend)
 
     Lx = 1.0
     Δx = Lx / nx
     Δy = Lx / ny
 
-    x = range(0, stop=Lx, length=nx)
-
     grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(Lx, Lx), dims=(nx, ny))
 
     # Courant number

examples/3D/3D_linear.jl

@@ -0,0 +1,74 @@
+using FiniteDiffWENO5
+using GLMakie
+
+function main(; nx=50, ny=50, nz=50)
+
+    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)
+
+    # Courant number
+    CFL = 0.7
+    period = 1
+
+    # 3D grid
+    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)
+
+    vx0 = ones(size(X))
+    vy0 = ones(size(Y))
+    vz0 = zeros(size(X)) # Rotation in XY plane only
+
+    v = (; x=vx0, y=vy0, z=vz0)
+
+    x0 = 1/4
+    c = 0.08
+
+    u0 = zeros(ny, nx, nz)
+    for I in CartesianIndices((ny, nx, nz))
+        u0[I] = exp(-((X[I]-x0)^2 + (Y[I]-x0)^2 + (Z[I]-x0)^2) / c^2)
+    end
+
+    u = copy(u0)
+    weno = WENOScheme(u; boundary=(2, 2, 2, 2, 2, 2), stag=false, multithreading=true)
+
+    Δt = CFL * min(Δx, Δy, Δz)^(5/3)
+    tmax = period * L / max(maximum(abs.(vx0)), maximum(abs.(vy0)), maximum(abs.(vz0)))
+    t = 0
+    counter = 0
+
+    f = Figure(size = (800, 600))
+    ax = Axis(f[1, 1], title = "t = $(round(t, digits=2))")
+
+    u_obser = Observable(u[:, :, div(nz, 2)])
+
+    heatmap!(ax, u_obser, colormap = :viridis)
+    Colorbar(f[1, 2], label = "u")
+    display(f)
+
+    while t < tmax
+        WENO_step!(u, v, weno, Δt, Δx, Δy, Δz)
+
+        t += Δt
+        if t + Δt > tmax
+            Δt = tmax - t
+        end
+
+        if counter % 50 == 0
+            u_obser[] = u[:, :, div(nz, 2)]
+            ax.title = "t = $(round(t, digits=2))"
+            sleep(0.01)
+        end
+
+        counter += 1
+    end
+end
+
+main()
+

examples/3D/3D_linear_chmy.jl

@@ -0,0 +1,86 @@
+using FiniteDiffWENO5
+using Chmy
+using KernelAbstractions
+using GLMakie
+
+function main(;backend=CPU(), nx=50, ny=50, nz=50)
+
+    arch = Arch(backend)
+
+    Lx = 1.0
+    Δx = Lx / nx
+    Δy = Lx / ny
+    Δz = Lx / nz
+
+    grid = UniformGrid(arch; origin=(0.0, 0.0, 0.0), extent=(Lx, Lx, Lx), dims=(nx, ny, nz))
+
+    # Courant number
+    CFL = 0.7
+    period = 1
+
+    # 3D grid
+    x = range(0, length=nx, stop=Lx)
+    y = range(0, length=ny, stop=Lx)
+    z = range(0, length=nz, stop=Lx)
+
+    X = reshape(x, nx, 1, 1)
+    Y = reshape(y, 1, ny, 1)
+    Z = reshape(z, 1, 1, nz)
+
+    X3D = X .+ 0 .* Y .+ 0 .* Z
+    Y3D = 0 .* X .+ Y .+ 0 .* Z
+    Z3D = 0 .* X .+ 0 .* Y .+ Z
+
+    vx0 = ones(size(X3D))
+    vy0 = ones(size(Y3D))
+    vz0 = zeros(size(Z3D)) # Rotation in XY plane only
+
+    v = (; x=vx0, y=vy0, z=vz0)
+
+    x0 = 1/4
+    c = 0.08
+
+    u0 = zeros(ny, nx, nz)
+    for I in CartesianIndices((ny, nx, nz))
+        u0[I] = exp(-((X3D[I]-x0)^2 + (Y3D[I]-x0)^2 + (Z3D[I]-x0)^2) / c^2)
+    end
+
+    u = copy(u0)
+    weno = WENOScheme(u; boundary=(2, 2, 2, 2, 2, 2), stag=false, multithreading=true)
+
+    Δt = CFL * min(Δx, Δy, Δz)^(5/3)
+    tmax = period * Lx / max(maximum(abs.(vx0)), maximum(abs.(vy0)), maximum(abs.(vz0)))
+    t = 0
+    counter = 0
+
+    f = Figure(size = (800, 600))
+    ax = Axis(f[1, 1], title = "t = $(round(t, digits=2))")
+
+    u_obser = Observable(u[:, :, div(nz, 2)])
+
+    heatmap!(ax, u_obser, colormap = :viridis)
+    Colorbar(f[1, 2], label = "u")
+    display(f)
+
+    while t < tmax
+        WENO_step!(u, v, weno, Δt, Δx, Δy, Δz)
+
+        t += Δt
+        if t + Δt > tmax
+            Δt = tmax - t
+        end
+
+        if counter % 10 == 0
+            u_obser[] = u[:, :, div(nz, 2)]
+            ax.title = "t = $(round(t, digits=2))"
+            sleep(0.01)
+        end
+
+        counter += 1
+    end
+
+    return u
+end
+
+u = main()
+

ext/ChmyExt1D.jl

@@ -112,11 +112,11 @@ module ChmyExt1D
     """
         WENO_step!(u::T_field, v::NamedTuple{names, <:Tuple{<:T_field}}, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, grid::StructuredGrid, arch) where T_field <: AbstractField{<:Real} where names
 
-    Advance the solution `u` by one time step using the 3rd-order Runge-Kutta method with WENO5 spatial discretization using Chmy.jl fields.
+    Advance the solution `u` by one time step using the 3rd-order Runge-Kutta method with WENO5 spatial discretization using Chmy.jl fields in 1D.
 
     # Arguments
     - `u::T_field`: The current solution field to be updated in place.
-    - `v::NamedTuple{names, <:Tuple{<:T_field}}`: The velocity field (can be staggered or not based on `weno.stag`).
+    - `v::NamedTuple{names, <:Tuple{<:T_field}}`: The velocity field (can be staggered or not based on `weno.stag`). Needs to be a NamedTuple with field `:x`.
     - `weno::WENOScheme`: The WENO scheme structure containing necessary parameters and fields.
     - `Δt`: The time step size.
     - `Δx`: The spatial grid size.

ext/ChmyExt2D.jl

@@ -140,13 +140,13 @@ module ChmyExt2D
     end
 
     """
-        WENO_step!(u::T_field, v::NamedTuple{names, <:Tuple{<:T_field}}, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, grid::StructuredGrid, arch) where T_field <: AbstractField{<:Real} where names
+        WENO_step!(u::T_field, v, weno::FiniteDiffWENO5.WENOScheme, Δt, Δx, grid::StructuredGrid, arch) where T_field <: AbstractField{<:Real} where names
 
-    Advance the solution `u` by one time step using the 3rd-order Runge-Kutta method with WENO5 spatial discretization using Chmy.jl fields.
+    Advance the solution `u` by one time step using the 3rd-order Runge-Kutta method with WENO5 spatial discretization using Chmy.jl fields in 2D.
 
     # Arguments
     - `u::T_field`: The current solution field to be updated in place.
-    - `v::NamedTuple{names, <:Tuple{<:T_field}}`: The velocity field (can be staggered or not based on `weno.stag`).
+    - `v::NamedTuple{names, <:Tuple{<:T_field}}`: The velocity field (can be staggered or not based on `weno.stag`). Needs to be a NamedTuple with fields `:x` and `:y`.
     - `weno::WENOScheme`: The WENO scheme structure containing necessary parameters and fields.
     - `Δt`: The time step size.
     - `Δx`: The spatial grid size.