Skip to content

Commit

Permalink
Adopt new format
Browse files Browse the repository at this point in the history
  • Loading branch information
@vavrines
vavrines committed Nov 5, 2024
1 parent 1760b9f commit dd4a6a5
Show file tree
Hide file tree
Showing 57 changed files with 806 additions and 960 deletions.
9 changes: 8 additions & 1 deletion .JuliaFormatter.toml
View file
Original file line number Diff line number Diff line change
@@ -1 +1,8 @@
format_docstrings = true
always_for_in = true
remove_extra_newlines = true
always_use_return = true
whitespace_in_kwargs = false
format_docstrings = true
separate_kwargs_with_semicolon = true
disallow_single_arg_nesting = true
format_markdown = false
6 changes: 3 additions & 3 deletions dev/check_phi.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ Vr, Vs = ∂vandermonde_matrix(N, pl[:, 1], pl[:, 2])

pf, wf = triface_quadrature(N)
ψf = zeros(3, N + 1, Np)
for i = 1:3
for i in 1:3
ψf[i, :, :] .= vandermonde_matrix(N, pf[i, :, 1], pf[i, :, 2])
end

Expand All @@ -111,14 +111,14 @@ V1 = hcat(V[1, :], V[3, :], V[5, :], V[2, :], V[4, :], V[6, :]) |> permutedims

# Vandermonde -> flux points
Vf = zeros(9, 6)
for i = 1:3, j = 1:3
for i in 1:3, j in 1:3
Vf[(i-1)*3+j, :] .= ψf[i, j, :]
end
@test Vf ≈ py_Vf

# coefficients
ϕ1 = zeros(6, 9)
for i = 1:3, j = 1:3
for i in 1:3, j in 1:3
ϕ1[:, (i-1)*3+j] .= ϕ[i, j, :]
end
ϕ2 = hcat(ϕ1[1, :], ϕ1[3, :], ϕ1[5, :], ϕ1[2, :], ϕ1[4, :], ϕ1[6, :]) |> permutedims
Expand Down
50 changes: 24 additions & 26 deletions dev/curve_euler.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,14 +47,14 @@ for i in axes(u0, 1), j in axes(u0, 2), k in axes(u0, 3), l in axes(u0, 4)
u0[i, j, k, l, :] .= prim_conserve(prim, ks.gas.γ)
end

n1 = [[0.0, 0.0] for i = 1:ps.nr+1, j = 1:ps.nθ]
for i = 1:ps.nr+1, j = 1:ps.nθ
n1 = [[0.0, 0.0] for i in 1:ps.nr+1, j in 1:ps.nθ]
for i in 1:ps.nr+1, j in 1:ps.nθ
angle = sum(ps.dθ[1, 1:j-1]) + 0.5 * ps.dθ[1, j]
n1[i, j] .= [cos(angle), sin(angle)]
end

n2 = [[0.0, 0.0] for i = 1:ps.nr, j = 1:ps.nθ+1]
for i = 1:ps.nr, j = 1:ps.nθ+1
n2 = [[0.0, 0.0] for i in 1:ps.nr, j in 1:ps.nθ+1]
for i in 1:ps.nr, j in 1:ps.nθ+1
angle = π / 2 + sum(ps.dθ[1, 1:j-1])
n2[i, j] .= [cos(angle), sin(angle)]
end
Expand All @@ -69,22 +69,22 @@ function dudt!(du, u, p, t)
nsp = size(u, 3)

f = zeros(nx, ny, nsp, nsp, 4, 2)
for i in axes(f, 1), j in axes(f, 2), k = 1:nsp, l = 1:nsp
for i in axes(f, 1), j in axes(f, 2), k in 1:nsp, l in 1:nsp
fg, gg = euler_flux(u[i, j, k, l, :], γ)
for m = 1:4
for m in 1:4
f[i, j, k, l, m, :] .= inv(J[i, j][k, l]) * [fg[m], gg[m]]
end
end

u_face = zeros(nx, ny, 4, nsp, 4)
f_face = zeros(nx, ny, 4, nsp, 4, 2)
for i in axes(u_face, 1), j in axes(u_face, 2), l = 1:nsp, m = 1:4
for i in axes(u_face, 1), j in axes(u_face, 2), l in 1:nsp, m in 1:4
u_face[i, j, 1, l, m] = dot(u[i, j, l, :, m], ll)
u_face[i, j, 2, l, m] = dot(u[i, j, :, l, m], lr)
u_face[i, j, 3, l, m] = dot(u[i, j, l, :, m], lr)
u_face[i, j, 4, l, m] = dot(u[i, j, :, l, m], ll)

for n = 1:2
for n in 1:2
f_face[i, j, 1, l, m, n] = dot(f[i, j, l, :, m, n], ll)
f_face[i, j, 2, l, m, n] = dot(f[i, j, :, l, m, n], lr)
f_face[i, j, 3, l, m, n] = dot(f[i, j, l, :, m, n], lr)
Expand All @@ -93,7 +93,7 @@ function dudt!(du, u, p, t)
end

fx_interaction = zeros(nx + 1, ny, nsp, 4)
for i = 2:nx, j = 1:ny, k = 1:nsp
for i in 2:nx, j in 1:ny, k in 1:nsp
#=fw = @view fx_interaction[i, j, k, :]
uL = local_frame(u_face[i-1, j, 2, k, :], n1[i, j][1], n1[i, j][2])
uR = local_frame(u_face[i, j, 4, k, :], n1[i, j][1], n1[i, j][2])
Expand All @@ -105,7 +105,7 @@ function dudt!(du, u, p, t)
40dt .* (u_face[i, j, 4, k, :] - u_face[i-1, j, 2, k, :])
end
fy_interaction = zeros(nx, ny + 1, nsp, 4)
for i = 1:nx, j = 2:ny, k = 1:nsp
for i in 1:nx, j in 2:ny, k in 1:nsp
#=fw = @view fy_interaction[i, j, k, :]
uL = local_frame(u_face[i, j-1, 3, k, :], n2[i, j][1], n2[i, j][2])
uR = local_frame(u_face[i, j, 1, k, :], n2[i, j][1], n2[i, j][2])
Expand All @@ -118,15 +118,15 @@ function dudt!(du, u, p, t)
end

rhs1 = zeros(nx, ny, nsp, nsp, 4)
for i = 1:nx, j = 1:ny, k = 1:nsp, l = 1:nsp, m = 1:4
for i in 1:nx, j in 1:ny, k in 1:nsp, l in 1:nsp, m in 1:4
rhs1[i, j, k, l, m] = dot(f[i, j, :, l, m, 1], lpdm[k, :])
end
rhs2 = zeros(nx, ny, nsp, nsp, 4)
for i = 1:nx, j = 1:ny, k = 1:nsp, l = 1:nsp, m = 1:4
for i in 1:nx, j in 1:ny, k in 1:nsp, l in 1:nsp, m in 1:4
rhs2[i, j, k, l, m] = dot(f[i, j, k, :, m, 2], lpdm[l, :])
end

for i = 2:nx-1, j = 2:ny-1, k = 1:nsp, l = 1:nsp, m = 1:4
for i in 2:nx-1, j in 2:ny-1, k in 1:nsp, l in 1:nsp, m in 1:4
#=fxL = (inv(ps.Ji[i, j][4, l]) * n1[i, j])[1] * fx_interaction[i, j, l, m]
fxR = (inv(ps.Ji[i, j][2, l]) * n1[i+1, j])[1] * fx_interaction[i+1, j, l, m]
fyL = (inv(ps.Ji[i, j][1, k]) * n2[i, j])[2] * fy_interaction[i, j, l, m]
Expand All @@ -140,14 +140,12 @@ function dudt!(du, u, p, t)
(fyR - f_face[i, j, 3, k, m, 2]) * dhr[l]
)=#

du[i, j, k, l, m] = -(
rhs1[i, j, k, l, m] +
rhs2[i, j, k, l, m] +
(fx_interaction[i, j, l, m] - f_face[i, j, 4, l, m, 1]) * dhl[k] +
(fx_interaction[i+1, j, l, m] - f_face[i, j, 2, l, m, 1]) * dhr[k] +
(fy_interaction[i, j, k, m] - f_face[i, j, 1, k, m, 2]) * dhl[l] +
(fy_interaction[i, j+1, k, m] - f_face[i, j, 3, k, m, 2]) * dhr[l]
)
du[i, j, k, l, m] = -(rhs1[i, j, k, l, m] +
rhs2[i, j, k, l, m] +
(fx_interaction[i, j, l, m] - f_face[i, j, 4, l, m, 1]) * dhl[k] +
(fx_interaction[i+1, j, l, m] - f_face[i, j, 2, l, m, 1]) * dhr[k] +
(fy_interaction[i, j, k, m] - f_face[i, j, 1, k, m, 2]) * dhl[l] +
(fy_interaction[i, j+1, k, m] - f_face[i, j, 3, k, m, 2]) * dhr[l])
end

return nothing
Expand All @@ -159,18 +157,18 @@ prob = ODEProblem(dudt!, u0, tspan, p)

dt = 0.002
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Midpoint(), save_everystep = false, adaptive = false, dt = dt)
itg = init(prob, Midpoint(); save_everystep=false, adaptive=false, dt=dt)

@showprogress for iter = 1:50
@showprogress for iter in 1:50
step!(itg)
end

contourf(ps.x, ps.y, itg.u[:, :, 2, 2, 1], aspect_ratio = 1, legend = true)
contourf(ps.x, ps.y, itg.u[:, :, 2, 2, 1]; aspect_ratio=1, legend=true)

sol = zeros(ps.nr, ps.nθ, 4)
for i = 1:ps.nr, j = 1:ps.nθ
for i in 1:ps.nr, j in 1:ps.nθ
sol[i, j, :] .= conserve_prim(itg.u[i, j, 2, 2, :], ks.gas.γ)
sol[i, j, 4] = 1 / sol[i, j, 4]
end

contourf(ps.x, ps.y, sol[:, :, 2], aspect_ratio = 1, legend = true)
contourf(ps.x, ps.y, sol[:, :, 2]; aspect_ratio=1, legend=true)
82 changes: 40 additions & 42 deletions dev/cylinder.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,24 +7,24 @@ pyplot()
cd(@__DIR__)

begin
set = Setup(
case = "cylinder",
space = "2d0f",
flux = "hll",
collision = "nothing",
interpOrder = 3,
limiter = "positivity",
boundary = "euler",
cfl = 0.1,
maxTime = 1.0,
set = Setup(;
case="cylinder",
space="2d0f",
flux="hll",
collision="nothing",
interpOrder=3,
limiter="positivity",
boundary="euler",
cfl=0.1,
maxTime=1.0,
)
ps = begin
ps0 = KitBase.CSpace2D(1.0, 6.0, 30, 0.0, π, 40, 0, 1)
deg = set.interpOrder - 1
FRPSpace2D(ps0, deg)
end
vs = nothing
gas = Gas(Kn = 1e-6, Ma = 1.2, K = 1.0)
gas = Gas(; Kn=1e-6, Ma=1.2, K=1.0)
ib = nothing

ks = SolverSet(set, ps0, vs, gas, ib)
Expand All @@ -36,14 +36,14 @@ for i in axes(u0, 1), j in axes(u0, 2), k in axes(u0, 3), l in axes(u0, 4)
u0[i, j, k, l, :] .= prim_conserve(prim, ks.gas.γ)
end

n1 = [[0.0, 0.0] for i = 1:ps.nr+1, j = 1:ps.nθ]
for i = 1:ps.nr+1, j = 1:ps.nθ
n1 = [[0.0, 0.0] for i in 1:ps.nr+1, j in 1:ps.nθ]
for i in 1:ps.nr+1, j in 1:ps.nθ
angle = sum(ps.dθ[1, 1:j-1]) + 0.5 * ps.dθ[1, j]
n1[i, j] .= [cos(angle), sin(angle)]
end

n2 = [[0.0, 0.0] for i = 1:ps.nr, j = 1:ps.nθ+1]
for i = 1:ps.nr, j = 1:ps.nθ+1
n2 = [[0.0, 0.0] for i in 1:ps.nr, j in 1:ps.nθ+1]
for i in 1:ps.nr, j in 1:ps.nθ+1
angle = π / 2 + sum(ps.dθ[1, 1:j-1])
n2[i, j] .= [cos(angle), sin(angle)]
end
Expand All @@ -58,22 +58,22 @@ function dudt!(du, u, p, t)
nsp = size(u, 3)

f = OffsetArray{Float64}(undef, 1:nx, 0:ny+1, nsp, nsp, 4, 2)
for i in axes(f, 1), j in axes(f, 2), k = 1:nsp, l = 1:nsp
for i in axes(f, 1), j in axes(f, 2), k in 1:nsp, l in 1:nsp
fg, gg = euler_flux(u[i, j, k, l, :], γ)
for m = 1:4
for m in 1:4
f[i, j, k, l, m, :] .= inv(J[i, j][k, l]) * [fg[m], gg[m]]
end
end

u_face = OffsetArray{Float64}(undef, 1:nx, 0:ny+1, 4, nsp, 4)
f_face = OffsetArray{Float64}(undef, 1:nx, 0:ny+1, 4, nsp, 4, 2)
for i in axes(u_face, 1), j in axes(u_face, 2), l = 1:nsp, m = 1:4
for i in axes(u_face, 1), j in axes(u_face, 2), l in 1:nsp, m in 1:4
u_face[i, j, 1, l, m] = dot(u[i, j, l, :, m], ll)
u_face[i, j, 2, l, m] = dot(u[i, j, :, l, m], lr)
u_face[i, j, 3, l, m] = dot(u[i, j, l, :, m], lr)
u_face[i, j, 4, l, m] = dot(u[i, j, :, l, m], ll)

for n = 1:2
for n in 1:2
f_face[i, j, 1, l, m, n] = dot(f[i, j, l, :, m, n], ll)
f_face[i, j, 2, l, m, n] = dot(f[i, j, :, l, m, n], lr)
f_face[i, j, 3, l, m, n] = dot(f[i, j, l, :, m, n], lr)
Expand All @@ -82,13 +82,13 @@ function dudt!(du, u, p, t)
end

fx_interaction = zeros(nx + 1, ny, nsp, 4)
for i = 2:nx, j = 1:ny, k = 1:nsp
for i in 2:nx, j in 1:ny, k in 1:nsp
fx_interaction[i, j, k, :] .=
0.5 .* (f_face[i-1, j, 2, k, :, 1] .+ f_face[i, j, 4, k, :, 1]) .-
dt .* (u_face[i, j, 4, k, :] - u_face[i-1, j, 2, k, :])
end

for j = 1:ny, k = 1:nsp
for j in 1:ny, k in 1:nsp
ul = local_frame(u_face[1, j, 4, k, :], n1[1, j][1], n1[1, j][2])
prim = conserve_prim(ul, γ)
pn = zeros(4)
Expand All @@ -103,7 +103,7 @@ function dudt!(du, u, p, t)

fg, gg = euler_flux(ub, γ)
fb = zeros(4)
for m = 1:4
for m in 1:4
fb[m] = (inv(ps.Ji[1, j][4, k])*[fg[m], gg[m]])[1]
end

Expand All @@ -112,30 +112,28 @@ function dudt!(du, u, p, t)
end

fy_interaction = zeros(nx, ny + 1, nsp, 4)
for i = 1:nx, j = 1:ny+1, k = 1:nsp
for i in 1:nx, j in 1:ny+1, k in 1:nsp
fy_interaction[i, j, k, :] .=
0.5 .* (f_face[i, j-1, 3, k, :, 2] .+ f_face[i, j, 1, k, :, 2]) .-
dt .* (u_face[i, j, 1, k, :] - u_face[i, j-1, 3, k, :])
end

rhs1 = zeros(nx, ny, nsp, nsp, 4)
for i = 1:nx, j = 1:ny, k = 1:nsp, l = 1:nsp, m = 1:4
for i in 1:nx, j in 1:ny, k in 1:nsp, l in 1:nsp, m in 1:4
rhs1[i, j, k, l, m] = dot(f[i, j, :, l, m, 1], lpdm[k, :])
end
rhs2 = zeros(nx, ny, nsp, nsp, 4)
for i = 1:nx, j = 1:ny, k = 1:nsp, l = 1:nsp, m = 1:4
for i in 1:nx, j in 1:ny, k in 1:nsp, l in 1:nsp, m in 1:4
rhs2[i, j, k, l, m] = dot(f[i, j, k, :, m, 2], lpdm[l, :])
end

for i = 1:nx-1, j = 1:ny, k = 1:nsp, l = 1:nsp, m = 1:4
du[i, j, k, l, m] = -(
rhs1[i, j, k, l, m] +
rhs2[i, j, k, l, m] +
(fx_interaction[i, j, l, m] - f_face[i, j, 4, l, m, 1]) * dhl[k] +
(fx_interaction[i+1, j, l, m] - f_face[i, j, 2, l, m, 1]) * dhr[k] +
(fy_interaction[i, j, k, m] - f_face[i, j, 1, k, m, 2]) * dhl[l] +
(fy_interaction[i, j+1, k, m] - f_face[i, j, 3, k, m, 2]) * dhr[l]
)
for i in 1:nx-1, j in 1:ny, k in 1:nsp, l in 1:nsp, m in 1:4
du[i, j, k, l, m] = -(rhs1[i, j, k, l, m] +
rhs2[i, j, k, l, m] +
(fx_interaction[i, j, l, m] - f_face[i, j, 4, l, m, 1]) * dhl[k] +
(fx_interaction[i+1, j, l, m] - f_face[i, j, 2, l, m, 1]) * dhr[k] +
(fy_interaction[i, j, k, m] - f_face[i, j, 1, k, m, 2]) * dhl[l] +
(fy_interaction[i, j+1, k, m] - f_face[i, j, 3, k, m, 2]) * dhr[l])
end

return nothing
Expand All @@ -147,10 +145,10 @@ prob = ODEProblem(dudt!, u0, tspan, p)

dt = 0.0002
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Midpoint(), save_everystep = false, adaptive = false, dt = dt)
itg = init(prob, Midpoint(); save_everystep=false, adaptive=false, dt=dt)

@showprogress for iter = 1:50#nt
for i = 1:ps.nr, k = 1:ps.deg+1, l = 1:ps.deg+1
@showprogress for iter in 1:50#nt
for i in 1:ps.nr, k in 1:ps.deg+1, l in 1:ps.deg+1
u1 = itg.u[i, 1, 4-k, 4-l, :]
ug1 = [u1[1], u1[2], -u1[3], u1[4]]
itg.u[i, 0, k, l, :] .= ug1
Expand All @@ -159,23 +157,23 @@ itg = init(prob, Midpoint(), save_everystep = false, adaptive = false, dt = dt)
ug2 = [u2[1], u2[2], -u2[3], u2[4]]
itg.u[i, ps.nθ+1, k, l, :] .= ug2
end
@inbounds for j = 1:ps.nθ÷2, k = 1:ps.deg+1, l = 1:ps.deg+1
@inbounds for j in 1:ps.nθ÷2, k in 1:ps.deg+1, l in 1:ps.deg+1
itg.u[ps.nr, j, k, l, :] .= itg.u[ps.nr-1, j, k, l, :]
end

step!(itg)
end

sol = zeros(ps.nr, ps.nθ, 4)
for i = 1:ps.nr, j = 1:ps.nθ
for i in 1:ps.nr, j in 1:ps.nθ
sol[i, j, :] .= conserve_prim(itg.u[i, j, 1, 1, :], ks.gas.γ)
sol[i, j, 4] = 1 / sol[i, j, 4]
end

contourf(
ps.x[1:ps.nr, 1:ps.nθ],
ps.y[1:ps.nr, 1:ps.nθ],
sol[:, :, 1],
aspect_ratio = 1,
legend = true,
sol[:, :, 1];
aspect_ratio=1,
legend=true,
)
Loading

0 comments on commit dd4a6a5

Please sign in to comment.