Skip to content

Commit 5da0850

Browse files
IddingsiteIddingsite
committed
Finish 2D
1 parent 526714e commit 5da0850

File tree

9 files changed

+575
-20
lines changed

9 files changed

+575
-20
lines changed

‎Project.toml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@ Chmy = "33a72cf0-4690-46d7-b987-06506c2248b9"
1111

1212
[extensions]
1313
ChmyExt1D = ["Chmy", "KernelAbstractions"]
14+
ChmyExt2D = ["Chmy", "KernelAbstractions"]
1415

1516
[compat]
1617
Chmy = "0.1"

‎examples/2D/2D_linear.jl‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -18,9 +18,8 @@ function main(;nx=400, ny=400)
1818
y = range(0, length=ny, stop= Lx)
1919
grid = (x .* ones(ny)', ones(nx) .* y')
2020

21-
w = π
22-
vx0 = w .* (grid[1] .- Lx/2)
23-
vy0 = -w .* (grid[2] .- Lx/2)
21+
vx0 = ones(nx, ny)
22+
vy0 = ones(nx, ny)
2423

2524
v = (;x=vy0, y=vx0)
2625

@@ -40,7 +39,7 @@ function main(;nx=400, ny=400)
4039
# grid size
4140
Δt = CFL*min(Δx, Δy)^(5/3)
4241

43-
tmax = period / (w/(2*π))
42+
tmax = period * Lx / max(maximum(abs.(vx0)), maximum(abs.(vy0)))
4443

4544
t = 0
4645
counter = 0

‎examples/2D/2D_linear_chmy.jl‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
using FiniteDiffWENO5
2+
using Chmy
3+
using KernelAbstractions
4+
using Plots
5+
6+
function main(;backend=CPU(), nx=400, ny=400, stag=true)
7+
8+
arch = Arch(backend)
9+
10+
Lx = 1.0
11+
Δx = Lx / nx
12+
Δy = Lx / ny
13+
14+
x = range(0, stop=Lx, length=nx)
15+
16+
grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(Lx, Lx), dims=(nx, ny))
17+
18+
# Courant number
19+
CFL = 0.7
20+
period = 1
21+
22+
# Grid x assumed defined
23+
x = range(0, length=nx, stop= Lx)
24+
y = range(0, length=ny, stop= Lx)
25+
grid_array = (x .* ones(ny)', ones(nx) .* y')
26+
27+
if stag
28+
vx0 = ones(nx+1, ny)
29+
vy0 = ones(nx, ny+1)
30+
else
31+
vx0 = ones(nx, ny)
32+
vy0 = ones(nx, ny)
33+
end
34+
35+
v = (;x=vy0, y=vx0)
36+
37+
x0 = 1/4
38+
c = 0.08
39+
40+
u0 = zeros(ny, nx)
41+
42+
for I in CartesianIndices((ny, nx))
43+
u0[I] = sign(exp(-((grid_array[1][I] - x0)^2 + (grid_array[2][I]' - x0)^2) / c^2) - 0.5) * 0.5 + 0.5
44+
end
45+
46+
u = Field(backend, grid, Center())
47+
set!(u, u0)
48+
49+
weno = WENOScheme(u; boundary=(1, 1, 1, 1), stag=stag, multithreading=true)
50+
51+
if stag
52+
v = VectorField(backend, grid)
53+
else
54+
v = (x=Field(arch, grid, Center()),
55+
y=Field(arch, grid, Center()))
56+
end
57+
58+
set!(v.x, vy0)
59+
set!(v.y, vx0)
60+
61+
# grid size
62+
Δt = CFL*min(Δx, Δy)^(5/3)
63+
64+
tmax = period * Lx / max(maximum(abs.(vx0)), maximum(abs.(vy0)))
65+
66+
t = 0
67+
counter = 0
68+
69+
mass_ini = sum(u0) * Δx * Δy
70+
71+
while t < tmax
72+
WENO_step!(u, v, weno, Δt, Δx, Δy, grid, arch)
73+
74+
75+
t += Δt
76+
77+
if t + Δt > tmax
78+
Δt = tmax - t
79+
end
80+
81+
if counter % 100 == 0
82+
83+
mass_ratio = (sum(u) * Δx * Δy) / mass_ini
84+
85+
heatmap(x, y, u, title="t = $(round(t, digits=2)) | mass = $mass_ratio", clims=(0,1))
86+
gui()
87+
end
88+
89+
counter += 1
90+
91+
end
92+
93+
end
94+
95+
main(backend=CPU(), nx=400, ny=400, stag=true)

‎examples/2D/2D_linear_rotation.jl‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
using FiniteDiffWENO5
2+
using Plots
3+
4+
function main(;nx=400, ny=400)
5+
6+
Lx = 1.0
7+
Δx = Lx / nx
8+
Δy = Lx / ny
9+
10+
x = range(0, stop=Lx, length=nx)
11+
12+
# Courant number
13+
CFL = 0.7
14+
period = 1
15+
16+
# Grid x assumed defined
17+
x = range(0, length=nx, stop= Lx)
18+
y = range(0, length=ny, stop= Lx)
19+
grid = (x .* ones(ny)', ones(nx) .* y')
20+
21+
w = π
22+
vx0 = w .* (grid[1] .- Lx/2)
23+
vy0 = -w .* (grid[2] .- Lx/2)
24+
25+
v = (;x=vy0, y=vx0)
26+
27+
x0 = 1/4
28+
c = 0.08
29+
30+
u0 = zeros(ny, nx)
31+
32+
for I in CartesianIndices((ny, nx))
33+
u0[I] = sign(exp(-((grid[1][I] - x0)^2 + (grid[2][I]' - x0)^2) / c^2) - 0.5) * 0.5 + 0.5
34+
end
35+
36+
u = copy(u0)
37+
weno = WENOScheme(u; boundary=(2, 2, 2, 2), stag=false, multithreading=true)
38+
39+
40+
# grid size
41+
Δt = CFL*min(Δx, Δy)^(5/3)
42+
43+
tmax = period / (w/(2*π))
44+
45+
t = 0
46+
counter = 0
47+
48+
while t < tmax
49+
WENO_step!(u, v, weno, Δt, Δx, Δy)
50+
51+
52+
t += Δt
53+
54+
if t + Δt > tmax
55+
Δt = tmax - t
56+
end
57+
58+
if counter % 100 == 0
59+
heatmap(x, y, u, title="t = $(round(t, digits=2))", clims=(0,1))
60+
gui()
61+
end
62+
63+
counter += 1
64+
65+
end
66+
67+
end
68+
69+
main(nx=400, ny=400)

‎examples/2D/2D_linear_rotation_chmy.jl‎

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
using FiniteDiffWENO5
2+
using Chmy
3+
using KernelAbstractions
4+
using Plots
5+
6+
function main(;backend=CPU(), nx=400, ny=400)
7+
# backend=CPU()
8+
# nx=400
9+
# ny=400
10+
11+
arch = Arch(backend)
12+
13+
Lx = 1.0
14+
Δx = Lx / nx
15+
Δy = Lx / ny
16+
17+
x = range(0, stop=Lx, length=nx)
18+
19+
grid = UniformGrid(arch; origin=(0.0, 0.0), extent=(Lx, Lx), dims=(nx, ny))
20+
21+
# Courant number
22+
CFL = 0.7
23+
period = 1
24+
25+
# Grid x assumed defined
26+
x = range(0, length=nx, stop= Lx)
27+
y = range(0, length=ny, stop= Lx)
28+
grid_array = (x .* ones(ny)', ones(nx) .* y')
29+
30+
w = π
31+
vx0 = w .* (grid_array[1] .- Lx/2)
32+
vy0 = -w .* (grid_array[2] .- Lx/2)
33+
34+
v = (;x=vy0, y=vx0)
35+
36+
x0 = 1/4
37+
c = 0.08
38+
39+
u0 = zeros(ny, nx)
40+
41+
for I in CartesianIndices((ny, nx))
42+
u0[I] = sign(exp(-((grid_array[1][I] - x0)^2 + (grid_array[2][I]' - x0)^2) / c^2) - 0.5) * 0.5 + 0.5
43+
end
44+
45+
u = Field(backend, grid, Center())
46+
set!(u, u0)
47+
48+
weno = WENOScheme(u; boundary=(2, 2, 2, 2), stag=false, multithreading=true)
49+
50+
v = (x=Field(arch, grid, Center()),
51+
y=Field(arch, grid, Center()))
52+
53+
set!(v.x, vy0)
54+
set!(v.y, vx0)
55+
56+
# grid size
57+
Δt = CFL*min(Δx, Δy)^(5/3)
58+
59+
tmax = period / (w/(2*π))
60+
61+
t = 0
62+
counter = 0
63+
64+
while t < tmax
65+
WENO_step!(u, v, weno, Δt, Δx, Δy, grid, arch)
66+
67+
68+
t += Δt
69+
70+
if t + Δt > tmax
71+
Δt = tmax - t
72+
end
73+
74+
if counter % 100 == 0
75+
heatmap(x, y, u, title="t = $(round(t, digits=2))", clims=(0,1))
76+
gui()
77+
end
78+
79+
counter += 1
80+
81+
end
82+
83+
end
84+
85+
main(backend=CPU(), nx=400, ny=400)

‎ext/ChmyExt1D.jl‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@ module ChmyExt1D
99

1010

1111
"""
12-
WENOScheme(u::AbstractField{T, N}, grid; boundary=(2, 2), stag=true, weno_z=false, multithreading=false) where {T, N}
12+
WENOScheme(u::AbstractField{T, N}, grid; boundary=(2, 2), stag=true, multithreading=false) where {T, N}
1313
1414
Create a WENO scheme structure for the given field `u` on the specified `grid` using Chmy.jl.
1515
@@ -18,9 +18,9 @@ module ChmyExt1D
1818
- `grid::StructuredGrid`: The computational grid.
1919
- `boundary::NTuple{2N, Int}`: A tuple specifying the boundary conditions for each dimension (0: homogeneous Neumann, 1: homogeneous Dirichlet, 2: periodic). Default is periodic (2).
2020
- `stag::Bool`: Whether the grid is staggered (velocities on cell faces) or not (velocities on cell centers).
21-
- `weno_z::Bool`: Whether to use the WENO-Z formulation (Borges et al. 2008) or not. Default is true.
21+
- `multithreading::Bool`: Whether to use multithreading (only for 2D and 3D). Default is false.
2222
"""
23-
function WENOScheme(u::AbstractField{T, N}, grid; boundary=(2, 2), stag=true, weno_z=true) where {T, N}
23+
function WENOScheme(u::AbstractField{T, N}, grid; boundary=(2, 2), stag=true) where {T, N}
2424

2525
# check that boundary conditions are correctly defined
2626
@assert length(boundary) == 2N "Boundary conditions must be a tuple of length $(2N) for $(N)D data."
@@ -42,7 +42,7 @@ module ChmyExt1D
4242
TFlux = typeof(fl)
4343
TArray = typeof(du)
4444

45-
return WENOScheme{T, TArray, TFlux, N_boundary}(stag=stag, boundary=boundary, weno_z=weno_z, multithreading=multithreading, fl=fl, fr=fr, du=du, ut=ut)
45+
return WENOScheme{T, TArray, TFlux, N_boundary}(stag=stag, boundary=boundary, multithreading=multithreading, fl=fl, fr=fr, du=du, ut=ut)
4646
end
4747

4848
@kernel function WENO_flux_chmy_1D(fl, fr, u, boundary, nx, χ, γ, ζ, ϵ, g::StructuredGrid, O)
@@ -88,8 +88,8 @@ module ChmyExt1D
8888
u5 = u[iee]
8989
u6 = u[ieee]
9090

91-
fl.x[i] = FiniteDiffWENO5.weno5_reconstruction_upwind(u1, u2, u3, u4, u5, χ, γ, ζ, ϵ)
92-
fr.x[i] = FiniteDiffWENO5.weno5_reconstruction_downwind(u2, u3, u4, u5, u6, χ, γ, ζ, ϵ)
91+
fl[i] = FiniteDiffWENO5.weno5_reconstruction_upwind(u1, u2, u3, u4, u5, χ, γ, ζ, ϵ)
92+
fr[i] = FiniteDiffWENO5.weno5_reconstruction_downwind(u2, u3, u4, u5, u6, χ, γ, ζ, ϵ)
9393
end
9494

9595
@kernel function WENO_semi_discretisation_weno5_chmy!(du, fl, fr, v, stag, Δx_, g::StructuredGrid, O)
@@ -135,17 +135,17 @@ module ChmyExt1D
135135

136136
@unpack ut, du, fl, fr, stag, boundary, χ, γ, ζ, ϵ = weno
137137

138-
launch(arch, grid, WENO_flux_chmy_1D => (fl, fr, u, boundary, nx, χ, γ, ζ, ϵ, grid))
138+
launch(arch, grid, WENO_flux_chmy_1D => (fl.x, fr.x, u, boundary, nx, χ, γ, ζ, ϵ, grid))
139139
launch(arch, grid, WENO_semi_discretisation_weno5_chmy! => (du, fl, fr, v, stag, Δx_, grid))
140140

141141
ut .= @muladd u .- Δt .* du
142142

143-
launch(arch, grid, WENO_flux_chmy_1D => (fl, fr, ut, boundary, nx, χ, γ, ζ, ϵ, grid))
143+
launch(arch, grid, WENO_flux_chmy_1D => (fl.x, fr.x, ut, boundary, nx, χ, γ, ζ, ϵ, grid))
144144
launch(arch, grid, WENO_semi_discretisation_weno5_chmy! => (du, fl, fr, v, stag, Δx_, grid))
145145

146146
ut .= @muladd 0.75 .* u .+ 0.25 .* ut .- 0.25 .* Δt .* du
147147

148-
launch(arch, grid, WENO_flux_chmy_1D => (fl, fr, ut, boundary, nx, χ, γ, ζ, ϵ, grid))
148+
launch(arch, grid, WENO_flux_chmy_1D => (fl.x, fr.x, ut, boundary, nx, χ, γ, ζ, ϵ, grid))
149149
launch(arch, grid, WENO_semi_discretisation_weno5_chmy! => (du, fl, fr, v, stag, Δx_, grid))
150150

151151
u .= @muladd inv(3.0) .* u .+ 2.0/3.0 .* ut .- 2.0/3.0 .* Δt .* du

0 commit comments

Comments
 (0)