diff --git a/dev/cylinder.jl b/dev/cylinder.jl index 0843bc94..ab76e283 100644 --- a/dev/cylinder.jl +++ b/dev/cylinder.jl @@ -19,17 +19,17 @@ begin ps = CSpace2D(1.0, 6.0, 30, 0.0, π, 50, 1, 1) vs = VSpace2D(-10.0, 10.0, 48, -10.0, 10.0, 48) gas = Gas(5e-2, 3.0, 1.0, 1.0) - + prim0 = [1.0, 0.0, 0.0, 1.0] prim1 = [1.0, gas.Ma * sound_speed(1.0, gas.γ), 0.0, 1.0] fw = (args...) -> prim_conserve(prim1, gas.γ) - ff = function(args...) + ff = function (args...) prim = conserve_prim(fw(args...), gas.γ) h = maxwellian(vs.u, vs.v, prim) b = h .* gas.K / 2 / prim[end] return h, b end - bc = function(x, y) + bc = function (x, y) if abs(x^2 + y^2 - 1) < 1e-3 return prim0 else @@ -48,14 +48,23 @@ dt = timestep(ks, ctr, 0.0) nt = ks.set.maxTime ÷ dt |> Int @showprogress for iter = 1:1000#nt evolve!(ks, ctr, a1face, a2face, dt; bc = [:maxwell, :extra, :mirror, :mirror]) - update!(ks, ctr, a1face, a2face, dt, zeros(4); bc = [:maxwell, :extra, :mirror, :mirror]) + update!( + ks, + ctr, + a1face, + a2face, + dt, + zeros(4); + bc = [:maxwell, :extra, :mirror, :mirror], + ) for j = ks.ps.nθ÷2+1:ks.ps.nθ ctr[ks.ps.nr+1, j].w .= ks.ib.fw(6, 0) ctr[ks.ps.nr+1, j].prim .= conserve_prim(ctr[ks.ps.nr+1, j].w, ks.gas.γ) ctr[ks.ps.nr+1, j].sw .= 0.0 ctr[ks.ps.nr+1, j].h .= maxwellian(ks.vs.u, ks.vs.v, ctr[ks.ps.nr+1, j].prim) - ctr[ks.ps.nr+1, j].b .= ctr[ks.ps.nr+1, j].h .* ks.gas.K ./ 2 ./ ctr[ks.ps.nr+1, j].prim[end] + ctr[ks.ps.nr+1, j].b .= + ctr[ks.ps.nr+1, j].h .* ks.gas.K ./ 2 ./ ctr[ks.ps.nr+1, j].prim[end] end global t += dt diff --git a/dev/reconstruct.jl b/dev/reconstruct.jl index c94130e3..ea82675a 100644 --- a/dev/reconstruct.jl +++ b/dev/reconstruct.jl @@ -57,7 +57,7 @@ nt = ks.set.maxTime ÷ dt |> Int @inbounds Threads.@threads for i in eachindex(face) vn = ks.vSpace.u .* face[i].n[1] .+ ks.vSpace.v .* face[i].n[2] vt = ks.vSpace.v .* face[i].n[1] .- ks.vSpace.u .* face[i].n[2] - + if !(-1 in ps.faceCells[i, :]) KitBase.flux_kfvs!( face[i].fw, @@ -94,7 +94,7 @@ nt = ks.set.maxTime ÷ dt |> Int dt, face[i].len, ) - + face[i].fw .= KitBase.global_frame(face[i].fw, face[i].n[1], face[i].n[2]) end end @@ -140,7 +140,7 @@ nt = ks.set.maxTime ÷ dt |> Int ) end end - + for i in eachindex(ps.cellType) if ps.cellType[i] == 3 ids = ps.cellNeighbors[i, :] diff --git a/example/advection.jl b/example/advection.jl index e0187340..a848f424 100644 --- a/example/advection.jl +++ b/example/advection.jl @@ -33,7 +33,7 @@ anim = @animate for iter = 1:nt evolve!(ks, ctr, face, dt) update!(ks, ctr, face, dt, 0.0) - plot(ks, ctr, ylims=[-1, 1]) + plot(ks, ctr, ylims = [-1, 1]) end gif(anim, "advection.gif", fps = 45) diff --git a/example/sod.jl b/example/sod.jl index cfb6ba63..12cf5b1c 100644 --- a/example/sod.jl +++ b/example/sod.jl @@ -16,7 +16,7 @@ set = Setup( ps = PSpace1D(0.0, 1.0, 100, 1) vs = VSpace1D(-5.0, 5.0, 72) gas = Gas(1e-4, 0.0, 1.0, 2.0) -fw, ff, bc = ib_sod(set, ps, vs, gas) +fw, ff, bc = config_ib(set, ps, vs, gas) ib = IB1F(fw, ff, bc) ks = SolverSet(set, ps, vs, gas, ib) @@ -32,4 +32,4 @@ for iter = 1:nt end plot(ks, ctr) -plot_line(ks, ctr) \ No newline at end of file +plot_line(ks, ctr) diff --git a/src/Boundary/bc_balance.jl b/src/Boundary/bc_balance.jl index 4b67cc79..7e65ce05 100644 --- a/src/Boundary/bc_balance.jl +++ b/src/Boundary/bc_balance.jl @@ -57,11 +57,7 @@ function bc_balance!( return nothing end -function bc_balance!( - ctr::T, - ctr0::T, - ctr1::T, -) where {T<:ControlVolume1D4F} +function bc_balance!(ctr::T, ctr0::T, ctr1::T) where {T<:ControlVolume1D4F} @. ctr.w = 2.0 * ctr0.w - ctr1.w @. ctr.prim = 2.0 * ctr0.prim - ctr1.prim @. ctr.h0 = 2.0 * ctr0.h0 - ctr1.h0 diff --git a/src/Boundary/bc_extra.jl b/src/Boundary/bc_extra.jl index 6404749d..a1b85956 100644 --- a/src/Boundary/bc_extra.jl +++ b/src/Boundary/bc_extra.jl @@ -1,7 +1,11 @@ +# ============================================================ +# Extrapolation Functions +# ============================================================ + function bc_extra!( ctr::AbstractVector{T}, ng = 1::Integer; - dirc + dirc, ) where {T<:AbstractControlVolume1D} if Symbol(dirc) in (:xl, :xL) @@ -16,13 +20,13 @@ function bc_extra!( end return nothing - + end function bc_extra!( ctr::AbstractMatrix{T}, ng = 1::Integer; - dirc + dirc, ) where {T<:AbstractControlVolume2D} if Symbol(dirc) in (:xl, :xL) @@ -46,5 +50,5 @@ function bc_extra!( end return nothing - + end diff --git a/src/Boundary/bc_mirror.jl b/src/Boundary/bc_mirror.jl index 819c82ae..1c2ee9e9 100644 --- a/src/Boundary/bc_mirror.jl +++ b/src/Boundary/bc_mirror.jl @@ -1,7 +1,11 @@ +# ============================================================ +# Mirroring Functions +# ============================================================ + function bc_mirror!( ctr::AbstractMatrix{T}, ng = 1::Integer; - dirc + dirc, ) where {T<:AbstractControlVolume2D} if Symbol(dirc) in (:xl, :xL) @@ -28,11 +32,7 @@ function bc_mirror!( end -function bc_mirror!( - ctr::ControlVolume2D, - ctr0::ControlVolume2D, - dirc -) +function bc_mirror!(ctr::ControlVolume2D, ctr0::ControlVolume2D, dirc) copy_ctr!(ctr, ctr0) if Symbol(dirc) == :x @@ -48,11 +48,7 @@ function bc_mirror!( return nothing end -function bc_mirror!( - ctr::ControlVolume2D1F, - ctr0::ControlVolume2D1F, - dirc -) +function bc_mirror!(ctr::ControlVolume2D1F, ctr0::ControlVolume2D1F, dirc) copy_ctr!(ctr, ctr0) nu = size(ctr.f, 1) @@ -79,11 +75,7 @@ function bc_mirror!( return nothing end -function bc_mirror!( - ctr::ControlVolume2D2F, - ctr0::ControlVolume2D2F, - dirc -) +function bc_mirror!(ctr::ControlVolume2D2F, ctr0::ControlVolume2D2F, dirc) copy_ctr!(ctr, ctr0) nu = size(ctr.h, 1) diff --git a/src/Boundary/bc_period.jl b/src/Boundary/bc_period.jl index 5c3eccaf..f9b69597 100644 --- a/src/Boundary/bc_period.jl +++ b/src/Boundary/bc_period.jl @@ -1,8 +1,8 @@ # ============================================================ -# Period Boundary Operations +# Period Functions # ============================================================ -function bc_period!(ctr::AbstractVector, ng = 1) +function bc_period!(ctr::AbstractVector{T}, ng = 1) where {T<:AbstractControlVolume} nx = length(ctr) - 2 * ng for i = 1:ng @@ -13,7 +13,7 @@ function bc_period!(ctr::AbstractVector, ng = 1) return nothing end -function bc_period!(ctr::AbstractMatrix, ng; dirc) +function bc_period!(ctr::AbstractMatrix{T}, ng; dirc) where {T<:AbstractControlVolume} if dirc == :x nx = size(ctr, 1) - 2 * ng for j in axes(ctr, 2), i = 1:ng diff --git a/src/Boundary/boundary.jl b/src/Boundary/boundary.jl index ad9c20ad..8aad3713 100644 --- a/src/Boundary/boundary.jl +++ b/src/Boundary/boundary.jl @@ -1,5 +1,5 @@ # ============================================================ -# Boundary Treatments +# Boundary Methods # ============================================================ include("bc_extra.jl") diff --git a/src/Config/brio-wu.txt b/src/Config/brio-wu.txt deleted file mode 100644 index 002f7595..00000000 --- a/src/Config/brio-wu.txt +++ /dev/null @@ -1,45 +0,0 @@ -# ------------------------------------------------------------ -# Configuration of Brio-Wu MHD shock tube -# ------------------------------------------------------------ - -# setup -case = brio-wu -space = 1d4f1v -nSpecies = 2 -interpOrder = 2 -limiter = vanleer -cfl = 0.5 -maxTime = 0.1 - -# phase space -x0 = 0 -x1 = 1 -nx = 100 -pMeshType = uniform -nxg = 0 - -# velocity space -u0 = -5.0 -u1 = 5.0 -nu = 28 -vMeshType = rectangle -nug = 0 - -# gas -knudsen = 0.00001 -mach = 0.0 -prandtl = 1.0 -inK = 0.0 - -mi = 1 -ni = 0.5 -#me = 0.04 -me = 0.0005446623 -ne = 0.5 -lD = 0.01 -rL = 0.003 - -# electromagnetic field -sol = 100 -echi = 1 -bnu = 1 diff --git a/src/Config/cfg_briowu.jl b/src/Config/cfg_briowu.jl index a466417b..02be1d22 100644 --- a/src/Config/cfg_briowu.jl +++ b/src/Config/cfg_briowu.jl @@ -1,5 +1,7 @@ """ - ib_briowu(gam, uspace::T, mi, me) where {T<:AbstractArray{<:AbstractFloat,2}} + ib_briowu(set, ps, vs, gas) + ib_briowu(gam, mi, me, uspace) + ib_briowu(gam, mi, me, uspace, vspace) Initialize Brio-Wu MHD shock tube @@ -52,28 +54,28 @@ function ib_briowu( BR[2] = -1.0 lorenzR = zeros(3, 2) - fw = function(x) + fw = function (x) if x <= (ps.x0 + ps.x1) / 2 return wL else return wR end end - fE = function(x) + fE = function (x) if x <= (ps.x0 + ps.x1) / 2 return EL else return ER end end - fB = function(x) + fB = function (x) if x <= (ps.x0 + ps.x1) / 2 return BL else return BR end end - fL = function(x) + fL = function (x) if x <= (ps.x0 + ps.x1) / 2 return lorenzL else @@ -81,7 +83,7 @@ function ib_briowu( end end - bc = function(x) + bc = function (x) if x <= (ps.x0 + ps.x1) / 2 return primL else @@ -112,7 +114,7 @@ function ib_briowu( (primR[3, j]^2 + primR[4, j]^2 + 2.0 / (2.0 * primR[end, j])) .* h0R[:, j] end - ff = function(x) + ff = function (x) if x <= (ps.x0 + ps.x1) / 2 return h0L, h1L, h2L, h3L else @@ -138,7 +140,7 @@ function ib_briowu( h2R[:, :, j] .= (primR[4, j]^2 + 1.0 / (2.0 * primR[end, j])) .* h0R[:, :, j] end - ff = function(x) + ff = function (x) if x <= (ps.x0 + ps.x1) / 2 return h0L, h1L, h2L else @@ -153,7 +155,7 @@ function ib_briowu( end -function ib_briowu(gam, mi, me, uspace::T) where {T<:AbstractArray{<:AbstractFloat,2}} +function ib_briowu(gam, mi, me, uspace::AbstractMatrix{T}) where {T<:AbstractFloat} # upstream primL = zeros(5, 2) diff --git a/src/Config/cfg_cavity.jl b/src/Config/cfg_cavity.jl index 248f64e5..e9515163 100644 --- a/src/Config/cfg_cavity.jl +++ b/src/Config/cfg_cavity.jl @@ -1,7 +1,5 @@ """ - 2d0f0v: ib_cavity(gam, Um, Vm, Tm) where {T<:AbstractArray{<:AbstractFloat,2}} - 2d1f2v: ib_cavity(gam, Um, Vm, Tm, u::T, v::T) where {T<:AbstractArray{<:AbstractFloat,2}} - 2d2f2v: ib_cavity(gam, Um, Vm, Tm, u::T, v::T, K) where {T<:AbstractArray{<:AbstractFloat,2}} + ib_cavity(set, ps, vs, gas, Um = 0.15, Vm = 0.0, Tm = 1.0) Initialize lid-driven cavity """ @@ -27,11 +25,11 @@ function ib_cavity( primL = deepcopy(prim) primR = deepcopy(prim) - fw = function(args...) + fw = function (args...) return w end - bc = function(x, y) + bc = function (x, y) if y == ps.y1 return primU else @@ -42,13 +40,13 @@ function ib_cavity( if set.space[1:4] == "2d0f" return fw, bc elseif set.space == "2d1f2v" - ff = function(args...) + ff = function (args...) return h end return fw, ff, bc elseif set.space == "2d2f2v" - ff = function(args...) + ff = function (args...) return h, b end diff --git a/src/Config/cfg_rh.jl b/src/Config/cfg_rh.jl index 5b328a7c..e330ae11 100644 --- a/src/Config/cfg_rh.jl +++ b/src/Config/cfg_rh.jl @@ -1,8 +1,6 @@ """ - 1d0f: ib_rh(MaL, gam) - 1d1f1v: ib_rh(MaL, gam, u::T) where {T<:AbstractArray{<:AbstractFloat,1}} - 1d2f1v: ib_rh(MaL, gam, u::T, K) where {T<:AbstractArray{<:AbstractFloat,1}} - 1d1f3v: ib_rh(MaL, gam, u::T, v::T, w::T) where {T<:AbstractArray{<:AbstractFloat,3}} + ib_rh(set, ps, vs, gas) + ib_rh(Ma, gam) Initialize Rankine-Hugoniot relation """ @@ -35,7 +33,7 @@ function ib_rh( wL = prim_conserve(primL, gam) wR = prim_conserve(primR, gam) - fw = function(x) + fw = function (x) if x <= (ps.x0 + ps.x1) / 2 return wL else @@ -43,7 +41,7 @@ function ib_rh( end end - bc = function(x) + bc = function (x) if x <= (ps.x0 + ps.x1) / 2 return primL else @@ -54,7 +52,7 @@ function ib_rh( if set.space[1:4] == "1d0f" return fw, bc elseif set.space == "1d1f1v" - ff = function(x) + ff = function (x) w = fw(x) prim = conserve_prim(w, gas.γ) h = maxwellian(vs.u, prim) @@ -63,7 +61,7 @@ function ib_rh( return fw, ff, bc elseif set.space == "1d2f1v" - ff = function(x) + ff = function (x) w = fw(x) prim = conserve_prim(w, gas.γ) h = maxwellian(vs.u, prim) @@ -73,7 +71,7 @@ function ib_rh( return fw, ff, bc elseif set.space == "1d1f3v" - ff = function(x) + ff = function (x) w = fw(x) prim = conserve_prim(w, gas.γ) h = maxwellian(vs.u, vs.v, vs.w, prim) @@ -99,7 +97,7 @@ function ib_rh( wL = mixture_prim_conserve(primL, gam) wR = mixture_prim_conserve(primR, gam) - fw = function(x) + fw = function (x) if x <= (ps.x0 + ps.x1) / 2 return wL else @@ -107,7 +105,7 @@ function ib_rh( end end - bc = function(x) + bc = function (x) if x <= (ps.x0 + ps.x1) / 2 return primL else @@ -125,7 +123,7 @@ function ib_rh( bR[:, j] .= hR[:, j] .* gas.K ./ (2.0 .* primR[end, j]) end - ff = function(x) + ff = function (x) if x <= (ps.x0 + ps.x1) / 2 return hL, bL else diff --git a/src/Config/cfg_sod.jl b/src/Config/cfg_sod.jl index c50cb864..3bc50cc4 100644 --- a/src/Config/cfg_sod.jl +++ b/src/Config/cfg_sod.jl @@ -1,8 +1,5 @@ """ - 1d0f0v: ib_sod(γ) - 1d1f1v: ib_sod(γ, u::T) where {T<:AbstractArray{<:AbstractFloat,1}} - 1d1f3v: ib_sod(γ, u::T, v::T, w::T) where {T<:AbstractArray{<:AbstractFloat,3}} - 1d2f1v: ib_sod(γ, u::T, K) where {T<:AbstractArray{<:AbstractFloat,1}} + ib_sod(set, ps, vs, gas) Initialize Sod shock tube """ @@ -21,7 +18,7 @@ function ib_sod( wL = prim_conserve(primL, gas.γ) wR = prim_conserve(primR, gas.γ) - fw = function(x) + fw = function (x) if x <= (ps.x0 + ps.x1) / 2 return wL else @@ -29,7 +26,7 @@ function ib_sod( end end - bc = function(x) + bc = function (x) if x <= (ps.x0 + ps.x1) / 2 return primL else @@ -40,7 +37,7 @@ function ib_sod( if set.space[1:4] == "1d0f" return fw, bc elseif set.space == "1d1f1v" - ff = function(x) + ff = function (x) w = fw(x) prim = conserve_prim(w, gas.γ) h = maxwellian(vs.u, prim) @@ -49,7 +46,7 @@ function ib_sod( return fw, ff, bc elseif set.space == "1d2f1v" - ff = function(x) + ff = function (x) w = fw(x) prim = conserve_prim(w, gas.γ) h = maxwellian(vs.u, prim) @@ -59,7 +56,7 @@ function ib_sod( return fw, ff, bc elseif set.space == "1d1f3v" - ff = function(x) + ff = function (x) w = fw(x) prim = conserve_prim(w, gas.γ) h = maxwellian(vs.u, vs.v, vs.w, prim) @@ -82,7 +79,7 @@ function ib_sod( wL = mixture_prim_conserve(primL, gas.γ) wR = mixture_prim_conserve(primR, gas.γ) - fw = function(x) + fw = function (x) if x <= (ps.x0 + ps.x1) / 2 return wL else @@ -90,7 +87,7 @@ function ib_sod( end end - bc = function(x) + bc = function (x) if x <= (ps.x0 + ps.x1) / 2 return primL else @@ -104,7 +101,7 @@ function ib_sod( hL = mixture_maxwellian(vs.u, primL) hR = mixture_maxwellian(vs.u, primR) - ff = function(x) + ff = function (x) if x <= (ps.x0 + ps.x1) / 2 return hL else @@ -123,7 +120,7 @@ function ib_sod( bR[:, j] .= hR[:, j] .* K ./ (2.0 .* primR[end, j]) end - ff = function(x) + ff = function (x) if x <= (ps.x0 + ps.x1) / 2 return hL, bL else diff --git a/src/Config/config.jl b/src/Config/config.jl index 08ed9689..26932139 100644 --- a/src/Config/config.jl +++ b/src/Config/config.jl @@ -2,7 +2,8 @@ # Initial & Boundary Conditions of Specific Problems # ============================================================ -export ib_rh, +export config_ib, + ib_rh, ib_sod, ib_briowu, ib_cavity @@ -11,3 +12,17 @@ include("cfg_rh.jl") include("cfg_sod.jl") include("cfg_briowu.jl") include("cfg_cavity.jl") + +function config_ib(args...; case = args[1].case) + func = begin + if case in ("shock", :shock) + eval(Symbol("ib_" * "rh")) + elseif case in ("brio-wu", Symbol("brio-wu")) + eval(Symbol("ib_" * "briowu")) + else + eval(Symbol("ib_" * string(case))) + end + end + + return func(args...) +end diff --git a/src/Config/shock.txt b/src/Config/shock.txt deleted file mode 100644 index 9e944699..00000000 --- a/src/Config/shock.txt +++ /dev/null @@ -1,35 +0,0 @@ -# ------------------------------------------------------------ -# Configuration of 1D shock structure -# ------------------------------------------------------------ - -# setup -case = shock -space = 1d2f1v -nSpecies = 1 -interpOrder = 1 -limiter = vanleer -cfl = 0.95 -maxTime = 250.0 - -# phase space -x0 = -35.0 -x1 = 35.0 -nx = 100 -pMeshType = uniform -nxg = 0 - -# velocity space -umin = -12.0 -umax = 12.0 -nu = 64 -vMeshType = rectangle -nug = 0 - -# gas -knudsen = 1.0 -mach = 2.0 -prandtl = 1.0 -inK = 0.0 -omega = 0.72 -alphaRef = 1.0 -omegaRef = 0.5 diff --git a/src/Data/data_array.jl b/src/Data/data_array.jl index d0038635..0b717975 100644 --- a/src/Data/data_array.jl +++ b/src/Data/data_array.jl @@ -1,3 +1,8 @@ +""" + slope_array(w; reduction = true) + +Generate array to store spatial slopes of solutions +""" slope_array(w::Number; kwargs...) = deepcopy(w) function slope_array(w::AbstractArray; reduction = true) diff --git a/src/Data/data_conversion.jl b/src/Data/data_conversion.jl index e7e093d8..b632e89e 100644 --- a/src/Data/data_conversion.jl +++ b/src/Data/data_conversion.jl @@ -1,7 +1,17 @@ +""" + symbolize(x) + +Generalized surrogate function of `Symbol` +""" symbolize(x::AbstractString) = Symbol(x) symbolize(x::AbstractArray{T}) where {T<:AbstractString} = [symbolize(y) for y in x] +""" + static_array(x) + +Transform to static array +""" function static_array(x::AbstractVector) y = MVector{length(x)}(collect(x)) @@ -30,7 +40,7 @@ function static_array(x::AbstractMatrix) return y end -function static_array(x::AbstractArray{<:Number,3}) +function static_array(x::AbstractArray{T,3}) where T y = MArray{Tuple{size(x, 1),size(x, 2),size(x, 3)}}(collect(x)) if x isa OffsetArray @@ -47,7 +57,7 @@ function static_array(x::AbstractArray{<:Number,3}) return y end -function static_array(x::AbstractArray{<:Number,4}) +function static_array(x::AbstractArray{T,4}) where T y = MArray{Tuple{size(x, 1),size(x, 2),size(x, 3),size(x, 4)}}(collect(x)) if x isa OffsetArray @@ -67,5 +77,10 @@ function static_array(x::AbstractArray{<:Number,4}) end -dynamic_array(x::AbstractArray) = x -dynamic_array(x::Number) = x +""" + dynamic_array(x) + +Transform to dynamic array +""" +dynamic_array(x) = x +dynamic_array(x::StaticArray) = Array(x) diff --git a/src/Flux/flux_boundary.jl b/src/Flux/flux_boundary.jl index 50ff3c34..f6686428 100644 --- a/src/Flux/flux_boundary.jl +++ b/src/Flux/flux_boundary.jl @@ -1,65 +1,30 @@ """ - flux_boundary_maxwell!( - fw::T1, - fh::T2, - fb::T2, - bc::T3, - h::T4, - b::T4, - u::T5, - ω::T5, - inK, - dt, - rot = 1, - ) where { - T1<:AbstractArray{<:AbstractFloat,1}, - T2<:AbstractArray{<:AbstractFloat,1}, - T3<:Array{<:Real,1}, - T4<:AbstractArray{<:AbstractFloat,1}, - T5<:AbstractArray{<:AbstractFloat,1}, - } - - flux_boundary_maxwell!( - fw::T1, - fh::T2, - fb::T2, - bc::T3, - h::T4, - b::T4, - u::T5, - v::T5, - ω::T5, - inK, - dt, - len, - rot = 1, - ) where { - T1<:AbstractArray{<:AbstractFloat,1}, - T2<:AbstractArray{<:AbstractFloat,2}, - T3<:Array{<:Real,1}, - T4<:AbstractArray{<:AbstractFloat,2}, - T5<:AbstractArray{<:AbstractFloat,2}, - } + flux_boundary_maxwell!(fw, bc, w, inK, γ, dt, rot) + flux_boundary_maxwell!(fw, bc, w, inK, γ, dt, len, rot) + flux_boundary_maxwell!(fw, fh, fb, bc, h, b, u, ω, inK, dt, rot) + flux_boundary_maxwell!(fw, ff, bc, f, u, v, ω, dt, len, rot) + flux_boundary_maxwell!(fw, fh, fb, bc, h, b, u, v, ω, inK, dt, len, rot) + flux_boundary_maxwell!(fw, ff, bc, f, u, v, w, ω, dt, area, rot) Maxwell's diffusive boundary flux -- @args: particle distribution functions and their slopes at left/right sides of interface +- @args: distribution functions and their slopes at left/right sides of interface - @args: particle velocity quadrature points and weights - @args: time step """ function flux_boundary_maxwell!( - fw::T1, - bc::T3, - w::T4, + fw::AbstractVector{T1}, + bc::AbstractVector{T2}, + w::AbstractVector{T3}, inK, γ, dt, rot, # 1 / -1 ) where { - T1<:AbstractVector{<:AbstractFloat}, - T3<:AbstractVector{<:Real}, - T4<:AbstractVector{<:Real}, + T1<:AbstractFloat, + T2<:Real, + T3<:Real, } # 1D continuum @assert length(bc) == 3 @@ -82,10 +47,10 @@ end #--- 1D2F1V ---# function flux_boundary_maxwell!( - fw::T1, + fw::AbstractVector{T1}, fh::T2, fb::T2, - bc::T3, + bc::AbstractVector{T3}, h::T4, b::T4, u::T5, @@ -94,11 +59,11 @@ function flux_boundary_maxwell!( dt, rot = 1, ) where { - T1<:AbstractArray{<:AbstractFloat,1}, - T2<:AbstractArray{<:AbstractFloat,1}, - T3<:Array{<:Real,1}, - T4<:AbstractArray{<:AbstractFloat,1}, - T5<:AbstractArray{<:AbstractFloat,1}, + T1<:AbstractFloat, + T2<:AbstractVector{<:AbstractFloat}, + T3<:Real, + T4<:AbstractVector{<:AbstractFloat}, + T5<:AbstractVector{<:AbstractFloat}, } # 1D2F1V @assert length(bc) == 3 @@ -131,18 +96,18 @@ end #--- 2D Continuum ---# function flux_boundary_maxwell!( - fw::T1, - bc::T3, - w::T4, + fw::AbstractVector{T1}, + bc::AbstractVector{T2}, + w::AbstractVector{T3}, inK, γ, dt, len, rot, ) where { - T1<:AbstractVector{<:AbstractFloat}, - T3<:AbstractVector{<:Real}, - T4<:AbstractVector{<:Real}, + T1<:Real, + T2<:Real, + T3<:Real, } @assert length(bc) == 4 @@ -167,10 +132,10 @@ end # 2D1F2V # ------------------------------------------------------------ function flux_boundary_maxwell!( - fw::T1, - fh::T2, - bc::T3, - h::T4, + fw::AbstractVector{T1}, + fh::AbstractMatrix{T2}, + bc::AbstractVector{T3}, + h::AbstractMatrix{T4}, u::T5, v::T5, ω::T5, @@ -178,11 +143,11 @@ function flux_boundary_maxwell!( len, rot = 1, ) where { - T1<:AbstractArray{<:AbstractFloat,1}, - T2<:AbstractArray{<:AbstractFloat,2}, - T3<:Array{<:Real,1}, - T4<:AbstractArray{<:AbstractFloat,2}, - T5<:AbstractArray{<:AbstractFloat,2}, + T1<:AbstractFloat, + T2<:AbstractFloat, + T3<:Real, + T4<:AbstractFloat, + T5<:AbstractMatrix{<:AbstractFloat}, } @assert length(bc) == 4 @@ -214,10 +179,10 @@ end # 2D2F2V # ------------------------------------------------------------ function flux_boundary_maxwell!( - fw::T1, + fw::AbstractVector{T1}, fh::T2, fb::T2, - bc::T3, + bc::AbstractVector{T3}, h::T4, b::T4, u::T5, @@ -228,11 +193,11 @@ function flux_boundary_maxwell!( len, rot = 1, ) where { - T1<:AbstractArray{<:AbstractFloat,1}, - T2<:AbstractArray{<:AbstractFloat,2}, - T3<:Array{<:Real,1}, - T4<:AbstractArray{<:AbstractFloat,2}, - T5<:AbstractArray{<:AbstractFloat,2}, + T1<:AbstractFloat, + T2<:AbstractMatrix{<:AbstractFloat}, + T3<:Real, + T4<:AbstractMatrix{<:AbstractFloat}, + T5<:AbstractMatrix{<:AbstractFloat}, } @assert length(bc) == 4 @@ -273,10 +238,10 @@ end # 1F3V # ------------------------------------------------------------ function flux_boundary_maxwell!( - fw::T1, - ff::T2, - bc::T3, - f::T4, + fw::AbstractVector{T1}, + ff::AbstractArray{T2,3}, + bc::AbstractVector{T3}, + f::AbstractArray{T4,3}, u::T5, v::T5, w::T5, @@ -285,10 +250,10 @@ function flux_boundary_maxwell!( area, rot = 1, ) where { - T1<:AbstractArray{<:AbstractFloat,1}, - T2<:AbstractArray{<:AbstractFloat,3}, - T3<:Array{<:Real,1}, - T4<:AbstractArray{<:AbstractFloat,3}, + T1<:AbstractFloat, + T2<:AbstractFloat, + T3<:Real, + T4<:AbstractFloat, T5<:AbstractArray{<:AbstractFloat,3}, } @@ -323,55 +288,29 @@ end """ - flux_boundary_specular!( - fw::T1, - ff::T2, - f::T3, - u::T4, - ω::T4, - dt, - ) where { - T1<:AbstractArray{<:Real,1}, - T2<:AbstractArray{<:AbstractFloat,1}, - T3<:AbstractArray{<:AbstractFloat,1}, - T4<:AbstractArray{<:AbstractFloat,1}, - } - - flux_boundary_specular!( - fw::T1, - fh::T2, - fb::T2, - h::T3, - b::T3, - u::T4, - ω::T4, - dt, - ) where { - T1<:AbstractArray{<:Real,1}, - T2<:AbstractArray{<:AbstractFloat,1}, - T3<:AbstractArray{<:AbstractFloat,1}, - T4<:AbstractArray{<:AbstractFloat,1}, - } + flux_boundary_specular!(fw, ff, f, u, ω, dt) + flux_boundary_specular!(fw, fh, fb, h, b, u, ω, dt) Specular reflection boundary flux -- @args: particle distribution functions and their slopes at left/right sides of interface -- @args: particle velocity quadrature points and weights +- @args: fluxes of conservative variables and distribution functions +- @args: distribution functions +- @args: velocity quadrature points and weights - @args: time step """ function flux_boundary_specular!( - fw::T1, - ff::T2, - f::T3, + fw::AbstractVector{T1}, + ff::AbstractVector{T2}, + f::AbstractVector{T3}, u::T4, ω::T4, dt, ) where { - T1<:AbstractArray{<:Real,1}, - T2<:AbstractArray{<:AbstractFloat,1}, - T3<:AbstractArray{<:AbstractFloat,1}, - T4<:AbstractArray{<:AbstractFloat,1}, + T1<:Real, + T2<:AbstractFloat, + T3<:AbstractFloat, + T4<:AbstractVector{<:AbstractFloat}, } # 1D1F1V fWall = similar(f) @@ -393,7 +332,7 @@ end # 1D2F1V # ------------------------------------------------------------ function flux_boundary_specular!( - fw::T1, + fw::AbstractVector{T1}, fh::T2, fb::T2, h::T3, @@ -402,10 +341,10 @@ function flux_boundary_specular!( ω::T4, dt, ) where { - T1<:AbstractArray{<:Real,1}, - T2<:AbstractArray{<:AbstractFloat,1}, - T3<:AbstractArray{<:AbstractFloat,1}, - T4<:AbstractArray{<:AbstractFloat,1}, + T1<:Real, + T2<:AbstractVector{<:AbstractFloat}, + T3<:AbstractVector{<:AbstractFloat}, + T4<:AbstractVector{<:AbstractFloat}, } hWall = similar(h) diff --git a/src/IO/io.jl b/src/IO/io.jl index b1468ff4..8ac6c9d4 100644 --- a/src/IO/io.jl +++ b/src/IO/io.jl @@ -26,7 +26,7 @@ function read_cfg(filename::T) where {T<:AbstractString} D[:boundary] = begin if parse(Int, D[:space][1]) == 1 [D[:boundary], D[:boundary]] - elseif parse(Int, D[:space][1]) == 2 + elseif parse(Int, D[:space][1]) == 2 [D[:boundary], D[:boundary], D[:boundary], D[:boundary]] end end @@ -190,7 +190,7 @@ end KS::X, ctr::Y, ) where {X<:AbstractSolverSet,Y<:AbstractArray{<:AbstractControlVolume,1}} - + # solution pltx = KS.pSpace.x[1:KS.pSpace.nx] plty = zeros(KS.pSpace.nx, 3) @@ -210,7 +210,7 @@ end label := "ρ" pltx, plty[:, 1] end - + if ctr[1].w isa AbstractArray @series begin label := "u" @@ -226,6 +226,7 @@ end c = get(plotattributes, :linewidth, :auto) nothing + end diff --git a/src/Math/math_entropy.jl b/src/Math/math_entropy.jl index 5cd3cbbb..8d12aa2f 100644 --- a/src/Math/math_entropy.jl +++ b/src/Math/math_entropy.jl @@ -3,7 +3,7 @@ Maxwell Boltzmann entropy """ -maxwell_boltzmann(f::T) where {T<:Real} = f * log(f) - f +maxwell_boltzmann(f) = f * log(f) - f """ @@ -31,7 +31,7 @@ maxwell_boltzmann_dual_prime(f) = exp(f) """ - kinetic_entropy(α::AbstractArray, m::AbstractArray, weights::AbstractVector) + kinetic_entropy(α, m, weights) Reconstruct mathematical entropy from Legendre dual """ diff --git a/src/Math/math_general.jl b/src/Math/math_general.jl index 08df7621..b4ad820e 100644 --- a/src/Math/math_general.jl +++ b/src/Math/math_general.jl @@ -1,5 +1,5 @@ """ - linspace(start, stop, n::T) where {T<:Integer} + linspace(start, stop, n) Python linspace function """ @@ -8,7 +8,7 @@ linspace(start, stop, n::T) where {T<:Integer} = """ - heaviside(x::Real) + heaviside(x) Heaviside step function """ @@ -16,19 +16,19 @@ heaviside(x::T) where {T<:Real} = ifelse(x >= 0, one(x), zero(x)) """ - fortsign(x::Real, y::Real) + fortsign(x, y) Fortran sign function """ -fortsign(x::T, y) where {T<:Real} = abs(x) * sign(y) +fortsign(x, y) = abs(x) * sign(y) """ - mat_split(m::AbstractArray) + mat_split(m) Split matrix into row vectors """ -function mat_split(m::T) where {T<:AbstractArray{<:Number,2}} +function mat_split(m::AbstractMatrix{T}) where {T<:Number} if length(m[1, :]) == 2 nx = eltype(m).([1.0 0.0]) @@ -47,12 +47,12 @@ end """ - central_diff(y::AbstractArray{<:Any,1}, x::AbstractArray{<:Any,1}) - central_diff(y::AbstractArray{<:Any,1}, dx::Any) + central_diff(y, x) + central_diff(y, dx) Central difference """ -function central_diff(y::T, x::T) where {T<:AbstractVector} +function central_diff(y::T, x::T) where {T<:AbstractVector{<:Number}} dy = zeros(eltype(y), axes(y)) @@ -70,7 +70,7 @@ function central_diff(y::T, x::T) where {T<:AbstractVector} end -function central_diff(y::T, dx) where {T<:AbstractVector} +function central_diff(y::AbstractVector{T}, dx) where T x = ones(eltype(y), axes(y)) .* dx dy = central_diff(y, x) @@ -79,11 +79,12 @@ end """ - central_diff2(y::AbstractArray{<:Any,1}, x::AbstractArray{<:Any,1}) + central_diff2(y, x) + central_diff2(y, dx) Central difference for second-order derivatives """ -function central_diff2(y::T, x::T) where {T<:AbstractVector} +function central_diff2(y::T, x::T) where {T<:AbstractVector{<:Number}} dy = zeros(eltype(y), axes(y)) i0 = eachindex(y) |> first @@ -96,7 +97,7 @@ function central_diff2(y::T, x::T) where {T<:AbstractVector} return dy end -function central_diff2(y::T, dx) where {T<:AbstractVector} +function central_diff2(y::AbstractVector{T}, dx) where T x = ones(eltype(y), axes(y)) .* dx dy = central_diff2(y, x) @@ -105,12 +106,12 @@ end """ - central_diff!(dy::AbstractArray{<:Any,1}, y::AbstractArray{<:Any,1}, x::AbstractArray{<:Any,1}) - central_diff!(dy::AbstractArray{<:Any,1}, y::AbstractArray{<:Any,1}, dx::Any) + central_diff!(dy, y, x) + central_diff!(dy, y, dx) Central difference """ -function central_diff!(dy::T1, y::T2, x::T2) where {T1<:AbstractVector,T2<:AbstractVector} +function central_diff!(dy::AbstractVector{T1}, y::T2, x::T2) where {T1,T2<:AbstractVector{<:Number}} @assert axes(dy) == axes(y) == axes(x) @@ -126,19 +127,19 @@ function central_diff!(dy::T1, y::T2, x::T2) where {T1<:AbstractVector,T2<:Abstr end -function central_diff!(dy::T1, y::T2, dx) where {T1<:AbstractVector,T2<:AbstractVector} +function central_diff!(dy::AbstractVector{T1}, y::AbstractVector{T2}, dx) where {T1,T2} x = ones(eltype(y), axes(y)) .* dx central_diff!(dy, y, x) end """ - central_diff2!(dy::T1, y::T2, x::T2) where {T1<:AbstractVector,T2<:AbstractVector} - central_diff2!(dy::T1, y::T2, dx) where {T1<:AbstractVector,T2<:AbstractVector} + central_diff2!(dy, y, x) + central_diff2!(dy, y, dx) Central difference """ -function central_diff2!(dy::T1, y::T2, x::T2) where {T1<:AbstractVector,T2<:AbstractVector} +function central_diff2!(dy::AbstractVector{T1}, y::T2, x::T2) where {T1,T2<:AbstractVector{<:Number}} @assert axes(dy) == axes(y) == axes(x) @@ -153,28 +154,23 @@ function central_diff2!(dy::T1, y::T2, x::T2) where {T1<:AbstractVector,T2<:Abst end -function central_diff2!(dy::T1, y::T2, dx) where {T1<:AbstractVector,T2<:AbstractVector} +function central_diff2!(dy::AbstractVector{T1}, y::AbstractVector{T2}, dx) where {T1,T2} x = ones(eltype(y), axes(y)) .* dx central_diff2!(dy, y, x) end """ - upwind_diff( - y::AbstractArray{<:Any,1}, - x::AbstractArray{<:Any,1}; - stream = :right::Symbol, - ) - - upwind_diff(y::AbstractArray{<:Any,1}, dx::Any; stream = :right::Symbol) + upwind_diff(y, x; stream) + upwind_diff(y, dx; stream) Upwind difference """ function upwind_diff( - y::AbstractArray{<:Any,1}, - x::AbstractArray{<:Any,1}; + y::AbstractVector{T1}, + x::AbstractVector{T2}; stream = :right::Symbol, -) +) where {T1,T2} dy = zeros(eltype(y), axes(y)) @@ -200,7 +196,7 @@ function upwind_diff( end -function upwind_diff(y::AbstractArray{<:Any,1}, dx::Any; stream = :right::Symbol) +function upwind_diff(y::AbstractVector{T}, dx; stream = :right::Symbol) where T x = ones(eltype(y), axes(y)) .* dx dy = upwind_diff(y, x, stream = stream) @@ -209,28 +205,17 @@ end """ - upwind_diff!( - dy::AbstractArray{<:Any,1}, - y::AbstractArray{<:Any,1}, - x::AbstractArray{<:Any,1}; - stream = :right::Symbol, - ) - - upwind_diff!( - dy::AbstractArray{<:Any,1}, - y::AbstractArray{<:Any,1}, - dx::Any; - stream = :right::Symbol, - ) + upwind_diff!(dy, y, x; stream) + upwind_diff!(dy, y, dx; stream) Upwind difference """ function upwind_diff!( - dy::AbstractArray{<:Any,1}, - y::AbstractArray{<:Any,1}, - x::AbstractArray{<:Any,1}; + dy::AbstractVector{T1}, + y::AbstractVector{T2}, + x::AbstractVector{T3}; stream = :right::Symbol, -) +) where {T1,T2,T3} @assert axes(dy) == axes(y) == axes(x) @@ -257,28 +242,31 @@ function upwind_diff!( end function upwind_diff!( - dy::AbstractArray{<:Any,1}, - y::AbstractArray{<:Any,1}, - dx::Any; + dy::AbstractVector{T1}, + y::AbstractVector{T2}, + dx; stream = :right::Symbol, -) +) where {T1,T2} + x = ones(eltype(y), axes(y)) .* dx upwind_diff!(dy, y, x, stream = stream) + end """ - unstruct_diff(u::AbstractArray{<:Any,1}, x::AbstractArray{<:Any,1}, nx::Integer; mode = :central::Symbol) - unstruct_diff(u::Function, x::AbstractArray{<:Any,2}, nx::Integer, dim::Integer; mode = :central::Symbol) + unstruct_diff(u, x, nx; mode) + unstruct_diff(u, x, nx, dim; mode) Finite difference for pseudo-unstructured mesh """ function unstruct_diff( - u::AbstractArray{<:Any,1}, - x::AbstractArray{<:Any,1}, + u::AbstractVector{T1}, + x::AbstractVector{T2}, nx::Integer; mode = :central::Symbol, -) +) where {T1,T2} + uu = reshape(u, (nx, :)) xx = reshape(x, (nx, :)) @@ -294,15 +282,17 @@ function unstruct_diff( end return reshape(dux, (1, :)) + end function unstruct_diff( u::Function, - x::AbstractArray{<:Any,1}, + x::AbstractVector{T}, nx::Integer, dim::Integer; mode = :central::Symbol, -) +) where T + uu = reshape(u.(x), (nx, :)) xx = reshape(x, (nx, :)) dux = zeros(eltype(x), axes(xx)) @@ -326,11 +316,12 @@ function unstruct_diff( end return dux[:] + end """ - lgwt(N::Integer, a::Real, b::Real) + lgwt(, a, b) Gauss Legendre integral for fast spectral method @@ -381,11 +372,11 @@ end """ - extract_last(a::AbstractArray, idx::Integer; mode=:view::Symbol) + extract_last(a, idx; mode) Extract subarray except the last column """ -function extract_last(a::AbstractArray, idx::Integer; mode = :view::Symbol) +function extract_last(a::AbstractArray{T}, idx::Integer; mode = :view::Symbol) where T if mode == :copy if ndims(a) == 2 diff --git a/src/Math/sphere_harmonics.jl b/src/Math/sphere_harmonics.jl index 20ac0e7d..5a065fcb 100644 --- a/src/Math/sphere_harmonics.jl +++ b/src/Math/sphere_harmonics.jl @@ -7,6 +7,7 @@ function ylmKCoefficient(l::Integer, m::Integer) return sqrt((2.0 * l + 1) / k / π) end + function ylmCosSinPolynomial(m::Integer, x::Variable, y::Variable) sum = 0.0 * (x + y) for j = 0:m÷2 @@ -16,6 +17,7 @@ function ylmCosSinPolynomial(m::Integer, x::Variable, y::Variable) return sum end + function ylmSinSinPolynomial(m::Integer, x::Variable, y::Variable) sum = 0.0 * (x + y) for j = 0:(m-1)÷2 @@ -25,8 +27,9 @@ function ylmSinSinPolynomial(m::Integer, x::Variable, y::Variable) return sum end + """ - ylm(l::Integer, m::Integer, x::Variable, y::Variable, z::Variable) + ylm(l, m, x, y, z) Calculation of the spherical harmonic for a given order (l,m) in Cartesian coordinates @@ -57,10 +60,11 @@ function ylm(l::Integer, m::Integer, x::Variable, y::Variable, z::Variable) end end + """ - rlylm(l::Integer, m::Integer, x::Variable, y::Variable, z::Variable) + rlylm(l, m, x, y, z) -Calculate r^l * Ylm(x,y,z) +Calculate `r^l * Ylm(x,y,z)`` """ function rlylm(l::Integer, m::Integer, x::Variable, y::Variable, z::Variable) p = ylm(l, m, x, y, z) @@ -77,11 +81,11 @@ end """ - eval_spherharmonic(points::T, L) where {T<:AbstractArray{<:Real,2}} + eval_spherharmonic(points, L) Evaluate spherical harmonics basis at given quadrature points """ -function eval_spherharmonic(points::T, L) where {T<:AbstractArray{<:Real,2}} +function eval_spherharmonic(points::AbstractMatrix{T}, L) where {T<:Real} ne = (L + 1)^2 nq = size(points, 1) m = zeros(ne, nq) diff --git a/src/Math/sphere_monomials.jl b/src/Math/sphere_monomials.jl index ba70020b..030bb9f2 100644 --- a/src/Math/sphere_monomials.jl +++ b/src/Math/sphere_monomials.jl @@ -2,6 +2,7 @@ degree_size(currDegree, spatialDim) = factorial(currDegree + spatialDim - 1) ÷ (factorial(currDegree) * factorial(spatialDim - 1)) + function basis_size(LMaxDegree, spatialDim) basisLen = 0 for idx_degree = 0:LMaxDegree @@ -10,6 +11,7 @@ function basis_size(LMaxDegree, spatialDim) return basisLen end + function power(basis, exponent) if exponent == 0 return 1.0 @@ -21,6 +23,7 @@ function power(basis, exponent) return result end + function monomial_basis(pointX, pointY, pointZ, polyDegree::T) where {T<:Integer} idx_vector = 1 spatialDim = 3 @@ -40,6 +43,7 @@ function monomial_basis(pointX, pointY, pointZ, polyDegree::T) where {T<:Integer return basisAtPt end + function monomial_basis(pointX, polyDegree::T) where {T<:Integer} idx_vector = 1 spatialDim = 1 @@ -53,6 +57,7 @@ function monomial_basis(pointX, polyDegree::T) where {T<:Integer} return basisAtPt end + function eval_sphermonomial(quadpts::AbstractVector, polyDegree::Integer) monomialBasis = zeros(basis_size(polyDegree, 1), size(quadpts, 1)) @@ -63,6 +68,7 @@ function eval_sphermonomial(quadpts::AbstractVector, polyDegree::Integer) return monomialBasis end + function eval_sphermonomial(quadpts::Matrix, polyDegree::Integer) monomialBasis = zeros(basis_size(polyDegree, 3), size(quadpts, 1)) diff --git a/src/Phase/creamer.jl b/src/Phase/creamer.jl index 2faa3d79..2d9a386e 100644 --- a/src/Phase/creamer.jl +++ b/src/Phase/creamer.jl @@ -3,11 +3,11 @@ # ============================================================ """ - slerp(pt1::T, pt2::T, n::I) where {T<:AbstractArray{<:Real,1},I<:Int} + slerp(pt1, pt2, n) Spherical linear interpolation """ -function slerp(pt1::T, pt2::T, n::I) where {T<:AbstractArray{<:Real,1},I<:Int} +function slerp(pt1::T, pt2::T, n::Integer) where {T<:AbstractArray{<:Real,1}} if norm(pt1 - pt2) < 1e-10 # same points return repeat(pt1, 1, n) # return n copies of that point end @@ -21,7 +21,7 @@ end """ - unique(Points::X, Triangles::Y) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:Int,2}} + unique(Points, Triangles) Cleaner for all duplicate (non unique) entries of quadrature points and triangles @@ -31,9 +31,9 @@ Cleaner for all duplicate (non unique) entries of quadrature points and triangle """ function unique( - Points::X, - Triangles::Y, -) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:Int,2}} + Points::AbstractMatrix{T1}, + Triangles::AbstractMatrix{T2}, +) where {T1<:Real,T2<:Integer} nPoints = size(Points)[2] nTriangles = size(Triangles)[2] @@ -126,7 +126,7 @@ function area( B::T, C::T, geometry = :plane::Symbol, -) where {T<:AbstractArray{<:Real,1}} +) where {T<:AbstractVector{<:Real}} if geometry == :plane alpha = angle(B, A, C) @@ -158,7 +158,7 @@ function angle( A::T, C::T, geometry = :plane::Symbol, -) where {T<:AbstractArray{<:Real,1}} +) where {T<:AbstractVector{<:Real}} if geometry == :plane u, v = A - B, C - A @@ -188,7 +188,7 @@ function distance( v1::T, v2::T, geometry = :plane::Symbol, -) where {T<:AbstractArray{<:Real,1}} +) where {T<:AbstractVector{<:Real}} if geometry == :plane return norm(v1 - v2) @@ -206,7 +206,7 @@ function distance( end -function muphi_xyz!(muphi::T, xyz::T) where {T<:AbstractArray{<:Real,2}} +function muphi_xyz!(muphi::T, xyz::T) where {T<:AbstractMatrix{<:Real}} n = size(xyz, 1) for i = 1:n xyz[i, 1] = sqrt(1 - muphi[i, 1]^2) * cos(muphi[i, 2]) @@ -216,7 +216,7 @@ function muphi_xyz!(muphi::T, xyz::T) where {T<:AbstractArray{<:Real,2}} end -function xyz_muphi!(xyz::T, muphi::T) where {T<:AbstractArray{<:Real,2}} +function xyz_muphi!(xyz::T, muphi::T) where {T<:AbstractMatrix{<:Real}} n = size(xyz, 1) for i = 1:n muphi[i, 1] = xyz[i, 3] diff --git a/src/Reconstruction/recon_filter.jl b/src/Reconstruction/recon_filter.jl index 2b276754..6b61a7e1 100644 --- a/src/Reconstruction/recon_filter.jl +++ b/src/Reconstruction/recon_filter.jl @@ -1,3 +1,16 @@ +# ------------------------------------------------------------ +# Filter functions +# ------------------------------------------------------------ + +""" + modal_filter!(u, args...; filter) + +Filter of modal solutions + +- @arg u: 1D modal solution +- @arg args...: filter parameters including strength, norm, etc. +- @arg filter: symbolic filter options including :l2, l2opt, :l1, :lasso, :exp, :houli +""" function modal_filter!( u::AbstractArray{T}, args...; @@ -31,7 +44,8 @@ function filter_l2!(u::AbstractMatrix{T}, args...) where {T<:AbstractFloat} λx, λξ = args[1:2] for j in axes(u, 2), i in axes(u, 1) - u[i, j] /= (1.0 + λx * (i - p0 + 1)^2 * (i - p0)^2 + λξ * (j - q0 + 1)^2 * (j - q0)^2) + u[i, j] /= + (1.0 + λx * (i - p0 + 1)^2 * (i - p0)^2 + λξ * (j - q0 + 1)^2 * (j - q0)^2) end return nothing @@ -71,7 +85,9 @@ function filter_l2opt!(u::AbstractMatrix{T}, args...) where {T<:AbstractFloat} η = η0 end - u[i, j] /= (1.0 + λx * (i - p0 + 1)^2 * (i - p0)^2 + λξ * (j - q0 + 1)^2 * (j - q0)^2 - η) + u[i, j] /= ( + 1.0 + λx * (i - p0 + 1)^2 * (i - p0)^2 + λξ * (j - q0 + 1)^2 * (j - q0)^2 - η + ) end end @@ -105,7 +121,7 @@ function filter_l1!(u::AbstractMatrix{T}, args...) where {T<:AbstractFloat} λ1, λ2 = args[1:2] ℓ = args[3] for j in axes(u, 2), i in axes(u, 1) - sc = 1.0 - (λ1 * i * (i - 1) + λ2 * j * (j - 1)) * ℓ[i, j] / (abs(u[i,j]) + 1e-8) + sc = 1.0 - (λ1 * i * (i - 1) + λ2 * j * (j - 1)) * ℓ[i, j] / (abs(u[i, j]) + 1e-8) if sc < 0.0 sc = 0.0 end diff --git a/src/Reconstruction/recon_limiter.jl b/src/Reconstruction/recon_limiter.jl index 289d22e2..f8589064 100644 --- a/src/Reconstruction/recon_limiter.jl +++ b/src/Reconstruction/recon_limiter.jl @@ -1,5 +1,5 @@ # ------------------------------------------------------------ -# Slope limiter functions +# Limiter functions # ------------------------------------------------------------ """ @@ -7,35 +7,41 @@ linear average of slopes """ -linear(sL::T, sR::T) where {T} = 0.5 * (sL + sR) +linear(sL, sR) = 0.5 * (sL + sR) """ vanleer(sL, sR) - vanleer(sL, s, sR) + vanleer(sL, s, sR, connect = 2) van Leer limiter """ -vanleer(sL::T, sR::T) where {T} = +vanleer(sL, sR) = (fortsign(1.0, sL) + fortsign(1.0, sR)) * abs(sL) * abs(sR) / (abs(sL) + abs(sR) + 1.e-7) #--- triangle ---# -function vanleer(sL::T, s::T, sR::T) where {T} - #=δ = [ - (fortsign(1.0, sL) + fortsign(1.0, s)) * abs(sL) * abs(s) / - (abs(sL) + abs(s) + 1.e-7), - (fortsign(1.0, s) + fortsign(1.0, sR)) * abs(s) * abs(sR) / - (abs(s) + abs(sR) + 1.e-7), - (fortsign(1.0, sL) + fortsign(1.0, sR)) * abs(sL) * abs(sR) / - (abs(sL) + abs(sR) + 1.e-7), - ]=# - δ = ( - (fortsign(1.0, sL) + fortsign(1.0, s)) * abs(sL) * abs(s) / - (abs(sL) + abs(s) + 1.e-7), - (fortsign(1.0, s) + fortsign(1.0, sR)) * abs(s) * abs(sR) / - (abs(s) + abs(sR) + 1.e-7), - ) +function vanleer(sL, s, sR, connect = 2) + δ = begin + if connect == 2 + ( + (fortsign(1.0, sL) + fortsign(1.0, s)) * abs(sL) * abs(s) / + (abs(sL) + abs(s) + 1.e-7), + (fortsign(1.0, s) + fortsign(1.0, sR)) * abs(s) * abs(sR) / + (abs(s) + abs(sR) + 1.e-7), + ) + elseif connect == 3 + ( + (fortsign(1.0, sL) + fortsign(1.0, s)) * abs(sL) * abs(s) / + (abs(sL) + abs(s) + 1.e-7), + (fortsign(1.0, s) + fortsign(1.0, sR)) * abs(s) * abs(sR) / + (abs(s) + abs(sR) + 1.e-7), + (fortsign(1.0, sL) + fortsign(1.0, sR)) * abs(sL) * abs(sR) / + (abs(sL) + abs(sR) + 1.e-7), + ) + end + end + id = findmin(abs.(δ))[2] return δ[id] @@ -44,24 +50,29 @@ end """ minmod(sL, sR) - minmod(sL, s, sR) + minmod(sL, s, sR, connect = 2) MinMod limiter """ -minmod(sL::T, sR::T) where {T} = - 0.5 * (fortsign(1.0, sL) + fortsign(1.0, sR)) * min(abs(sR), abs(sL)) +minmod(sL, sR) = 0.5 * (fortsign(1.0, sL) + fortsign(1.0, sR)) * min(abs(sR), abs(sL)) #--- triangle ---# -function minmod(sL::T, s::T, sR::T) where {T} - #=δ = [ - 0.5 * (fortsign(1.0, sL) + fortsign(1.0, s)) * min(abs(s), abs(sL)), - 0.5 * (fortsign(1.0, s) + fortsign(1.0, sR)) * min(abs(sR), abs(s)), - 0.5 * (fortsign(1.0, sL) + fortsign(1.0, sR)) * min(abs(sR), abs(sL)), - ]=# - δ = ( - 0.5 * (fortsign(1.0, sL) + fortsign(1.0, s)) * min(abs(s), abs(sL)), - 0.5 * (fortsign(1.0, s) + fortsign(1.0, sR)) * min(abs(sR), abs(s)), - ) +function minmod(sL, s, sR, connect = 2) + δ = begin + if connect == 2 + ( + 0.5 * (fortsign(1.0, sL) + fortsign(1.0, s)) * min(abs(s), abs(sL)), + 0.5 * (fortsign(1.0, s) + fortsign(1.0, sR)) * min(abs(sR), abs(s)), + ) + elseif connect == 3 + ( + 0.5 * (fortsign(1.0, sL) + fortsign(1.0, s)) * min(abs(s), abs(sL)), + 0.5 * (fortsign(1.0, s) + fortsign(1.0, sR)) * min(abs(sR), abs(s)), + 0.5 * (fortsign(1.0, sL) + fortsign(1.0, sR)) * min(abs(sR), abs(sL)), + ) + end + end + id = findmin(abs.(δ))[2] return δ[id] @@ -73,8 +84,7 @@ end SuperBee limiter """ -function superbee(sL::T, sR::T) where {T} - +function superbee(sL, sR) if sR >= 0.5 * sL && sR <= 2.0 * sL return 0.5 * (fortsign(1.0, sL) + fortsign(1.0, sR)) * max(abs(sL), abs(sR)) elseif sR < 0.5 * sL && sR > 2.0 * sL @@ -82,7 +92,6 @@ function superbee(sL::T, sR::T) where {T} else return 0.0 end - end @@ -91,23 +100,22 @@ end van Albaba limiter """ -vanalbaba(sL::T, sR::T) where {T} = (sL^2 * sR + sL * sR^2) / (sL^2 + sR^2 + 1.e-7) +vanalbaba(sL, sR) = (sL^2 * sR + sL * sR^2) / (sL^2 + sR^2 + 1.e-7) """ - weno5(wL2::T, wL1::T, wN::T, wR1::T, wR2::T) where {T} + weno5(wL2, wL1, wN, wR1, wR2) 5th-order WENO-JS interpolation """ -function weno5(wL2::T, wL1::T, wN::T, wR1::T, wR2::T) where {T} - +function weno5(wL2, wL1, wN, wR1, wR2) ϵ = 1e-6 β0 = 13.0 / 12.0 * (wN - 2.0 * wR1 + wR2)^2 + 1.0 / 4.0 * (3.0 * wN - 4.0 * wR1 + wR2)^2 β1 = 13.0 / 12.0 * (wL1 - 2.0 * wN + wR1)^2 + 1.0 / 4.0 * (wL1 - wR1)^2 β2 = 13.0 / 12.0 * (wL2 - 2.0 * wL1 + wN)^2 + 1.0 / 4.0 * (wL2 - 4.0 * wL1 + 3.0 * wN)^2 - #--- right interface ---# + # right interface dr0 = 0.3 dr1 = 0.6 dr2 = 0.1 @@ -126,7 +134,7 @@ function weno5(wL2::T, wL1::T, wN::T, wR1::T, wR2::T) where {T} wR = ωr0 * qr0 + ωr1 * qr1 + ωr2 * qr2 - #--- left interface ---# + # left interface dl0 = 0.1 dl1 = 0.6 dl2 = 0.3 @@ -146,5 +154,4 @@ function weno5(wL2::T, wL1::T, wN::T, wR1::T, wR2::T) where {T} wL = αl0 * ql0 + αl1 * ql1 + αl2 * ql2 return wL, wR - end diff --git a/src/Solver/evolve_boundary.jl b/src/Solver/evolve_boundary.jl index 975e7b0a..ae93a20f 100644 --- a/src/Solver/evolve_boundary.jl +++ b/src/Solver/evolve_boundary.jl @@ -24,11 +24,7 @@ function evolve_boundary!( vt = KS.vSpace.v .* a1face[1, j].n[1] .- KS.vSpace.u .* a1face[1, j].n[2] xc = (KS.ps.vertices[1, j, 1, 1] + KS.ps.vertices[1, j, 4, 1]) / 2 yc = (KS.ps.vertices[1, j, 1, 2] + KS.ps.vertices[1, j, 4, 2]) / 2 - bcL = local_frame( - KS.ib.bc(xc, yc), - a1face[1, j].n[1], - a1face[1, j].n[2], - ) + bcL = local_frame(KS.ib.bc(xc, yc), a1face[1, j].n[1], a1face[1, j].n[2]) flux_boundary_maxwell!( a1face[1, j].fw, a1face[1, j].fh, @@ -51,19 +47,11 @@ function evolve_boundary!( if bcs[2] == :maxwell @inbounds Threads.@threads for j = 1:ny - vn = - KS.vSpace.u .* a1face[nx+1, j].n[1] .+ - KS.vSpace.v .* a1face[nx+1, j].n[2] - vt = - KS.vSpace.v .* a1face[nx+1, j].n[1] .- - KS.vSpace.u .* a1face[nx+1, j].n[2] + vn = KS.vSpace.u .* a1face[nx+1, j].n[1] .+ KS.vSpace.v .* a1face[nx+1, j].n[2] + vt = KS.vSpace.v .* a1face[nx+1, j].n[1] .- KS.vSpace.u .* a1face[nx+1, j].n[2] xc = (KS.ps.vertices[nx, j, 2, 1] + KS.ps.vertices[nx, j, 3, 1]) / 2 yc = (KS.ps.vertices[nx, j, 2, 2] + KS.ps.vertices[nx, j, 3, 2]) / 2 - bcR = local_frame( - KS.ib.bc(xc, yc), - a1face[nx+1, j].n[1], - a1face[nx+1, j].n[2], - ) + bcR = local_frame(KS.ib.bc(xc, yc), a1face[nx+1, j].n[1], a1face[nx+1, j].n[2]) flux_boundary_maxwell!( a1face[nx+1, j].fw, a1face[nx+1, j].fh, @@ -79,11 +67,8 @@ function evolve_boundary!( dy[nx, j], -1, ) - a1face[nx+1, j].fw .= global_frame( - a1face[nx+1, j].fw, - a1face[nx+1, j].n[1], - a1face[nx+1, j].n[2], - ) + a1face[nx+1, j].fw .= + global_frame(a1face[nx+1, j].fw, a1face[nx+1, j].n[1], a1face[nx+1, j].n[2]) end end @@ -93,11 +78,7 @@ function evolve_boundary!( vt = KS.vSpace.v .* a2face[i, 1].n[1] .- KS.vSpace.u .* a2face[i, 1].n[2] xc = (KS.ps.vertices[i, 1, 1, 1] + KS.ps.vertices[i, 1, 2, 1]) / 2 yc = (KS.ps.vertices[i, 1, 1, 2] + KS.ps.vertices[i, 1, 2, 2]) / 2 - bcD = local_frame( - KS.ib.bc(xc, yc), - a2face[i, 1].n[1], - a2face[i, 1].n[2], - ) + bcD = local_frame(KS.ib.bc(xc, yc), a2face[i, 1].n[1], a2face[i, 1].n[2]) flux_boundary_maxwell!( a2face[i, 1].fw, @@ -121,19 +102,11 @@ function evolve_boundary!( if bcs[4] == :maxwell @inbounds Threads.@threads for i = 1:nx - vn = - KS.vSpace.u .* a2face[i, ny+1].n[1] .+ - KS.vSpace.v .* a2face[i, ny+1].n[2] - vt = - KS.vSpace.v .* a2face[i, ny+1].n[1] .- - KS.vSpace.u .* a2face[i, ny+1].n[2] + vn = KS.vSpace.u .* a2face[i, ny+1].n[1] .+ KS.vSpace.v .* a2face[i, ny+1].n[2] + vt = KS.vSpace.v .* a2face[i, ny+1].n[1] .- KS.vSpace.u .* a2face[i, ny+1].n[2] xc = (KS.ps.vertices[i, ny, 3, 1] + KS.ps.vertices[i, ny, 4, 1]) / 2 yc = (KS.ps.vertices[i, ny, 3, 2] + KS.ps.vertices[i, ny, 4, 2]) / 2 - bcU = local_frame( - KS.ib.bc(xc, yc), - a2face[i, ny+1].n[1], - a2face[i, ny+1].n[2], - ) + bcU = local_frame(KS.ib.bc(xc, yc), a2face[i, ny+1].n[1], a2face[i, ny+1].n[2]) flux_boundary_maxwell!( a2face[i, ny+1].fw, @@ -150,12 +123,9 @@ function evolve_boundary!( dx[i, ny], -1, ) - a2face[i, ny+1].fw .= global_frame( - a2face[i, ny+1].fw, - a2face[i, ny+1].n[1], - a2face[i, ny+1].n[2], - ) + a2face[i, ny+1].fw .= + global_frame(a2face[i, ny+1].fw, a2face[i, ny+1].n[1], a2face[i, ny+1].n[2]) end end -end \ No newline at end of file +end diff --git a/src/Solver/solver_evolution.jl b/src/Solver/solver_evolution.jl index f28d57a8..f52195e1 100644 --- a/src/Solver/solver_evolution.jl +++ b/src/Solver/solver_evolution.jl @@ -858,11 +858,7 @@ function evolve!( vt = KS.vSpace.v .* a1face[1, j].n[1] .- KS.vSpace.u .* a1face[1, j].n[2] xc = (KS.ps.vertices[1, j, 1, 1] + KS.ps.vertices[1, j, 4, 1]) / 2 yc = (KS.ps.vertices[1, j, 1, 2] + KS.ps.vertices[1, j, 4, 2]) / 2 - bcL = local_frame( - KS.ib.bc(xc, yc), - a1face[1, j].n[1], - a1face[1, j].n[2], - ) + bcL = local_frame(KS.ib.bc(xc, yc), a1face[1, j].n[1], a1face[1, j].n[2]) flux_boundary_maxwell!( a1face[1, j].fw, @@ -882,19 +878,11 @@ function evolve!( end if bcs[2] == :maxwell @inbounds Threads.@threads for j = 1:ny - vn = - KS.vSpace.u .* a1face[nx+1, j].n[1] .+ - KS.vSpace.v .* a1face[nx+1, j].n[2] - vt = - KS.vSpace.v .* a1face[nx+1, j].n[1] .- - KS.vSpace.u .* a1face[nx+1, j].n[2] + vn = KS.vSpace.u .* a1face[nx+1, j].n[1] .+ KS.vSpace.v .* a1face[nx+1, j].n[2] + vt = KS.vSpace.v .* a1face[nx+1, j].n[1] .- KS.vSpace.u .* a1face[nx+1, j].n[2] xc = (KS.ps.vertices[nx, j, 2, 1] + KS.ps.vertices[nx, j, 3, 1]) / 2 yc = (KS.ps.vertices[nx, j, 2, 2] + KS.ps.vertices[nx, j, 3, 2]) / 2 - bcR = local_frame( - KS.ib.bc(xc, yc), - a1face[nx+1, j].n[1], - a1face[nx+1, j].n[2], - ) + bcR = local_frame(KS.ib.bc(xc, yc), a1face[nx+1, j].n[1], a1face[nx+1, j].n[2]) flux_boundary_maxwell!( a1face[nx+1, j].fw, @@ -908,11 +896,8 @@ function evolve!( KS.ps.dy[nx, j], -1, ) - a1face[nx+1, j].fw .= global_frame( - a1face[nx+1, j].fw, - a1face[nx+1, j].n[1], - a1face[nx+1, j].n[2], - ) + a1face[nx+1, j].fw .= + global_frame(a1face[nx+1, j].fw, a1face[nx+1, j].n[1], a1face[nx+1, j].n[2]) end end if bcs[3] == :maxwell @@ -921,11 +906,7 @@ function evolve!( vt = KS.vSpace.v .* a2face[i, 1].n[1] .- KS.vSpace.u .* a2face[i, 1].n[2] xc = (KS.ps.vertices[i, 1, 1, 1] + KS.ps.vertices[i, 1, 2, 1]) / 2 yc = (KS.ps.vertices[i, 1, 1, 2] + KS.ps.vertices[i, 1, 2, 2]) / 2 - bcD = local_frame( - KS.ib.bc(xc, yc), - a2face[i, 1].n[1], - a2face[i, 1].n[2], - ) + bcD = local_frame(KS.ib.bc(xc, yc), a2face[i, 1].n[1], a2face[i, 1].n[2]) flux_boundary_maxwell!( a2face[i, 1].fw, @@ -945,19 +926,11 @@ function evolve!( end if bcs[4] == :maxwell @inbounds Threads.@threads for i = 1:nx - vn = - KS.vSpace.u .* a2face[i, ny+1].n[1] .+ - KS.vSpace.v .* a2face[i, ny+1].n[2] - vt = - KS.vSpace.v .* a2face[i, ny+1].n[1] .- - KS.vSpace.u .* a2face[i, ny+1].n[2] + vn = KS.vSpace.u .* a2face[i, ny+1].n[1] .+ KS.vSpace.v .* a2face[i, ny+1].n[2] + vt = KS.vSpace.v .* a2face[i, ny+1].n[1] .- KS.vSpace.u .* a2face[i, ny+1].n[2] xc = (KS.ps.vertices[i, ny, 3, 1] + KS.ps.vertices[i, ny, 4, 1]) / 2 yc = (KS.ps.vertices[i, ny, 3, 2] + KS.ps.vertices[i, ny, 4, 2]) / 2 - bcU = local_frame( - KS.ib.bc(xc, yc), - a2face[i, ny+1].n[1], - a2face[i, ny+1].n[2], - ) + bcU = local_frame(KS.ib.bc(xc, yc), a2face[i, ny+1].n[1], a2face[i, ny+1].n[2]) flux_boundary_maxwell!( a2face[i, ny+1].fw, @@ -971,11 +944,8 @@ function evolve!( KS.ps.dx[i, ny], -1, ) - a2face[i, ny+1].fw .= global_frame( - a2face[i, ny+1].fw, - a2face[i, ny+1].n[1], - a2face[i, ny+1].n[2], - ) + a2face[i, ny+1].fw .= + global_frame(a2face[i, ny+1].fw, a2face[i, ny+1].n[1], a2face[i, ny+1].n[2]) end end @@ -1285,7 +1255,11 @@ function evolve!( idx = ifelse(KS.ps.faceCells[i, 1] != -1, 1, 2) if KS.ps.cellType[KS.ps.faceCells[i, idx]] == 2 - bc = KitBase.local_frame(KS.ib.bc(KS.ps.faceCenter[i, 1], KS.ps.faceCenter[i, 2]), face[i].n[1], face[i].n[2]) + bc = KitBase.local_frame( + KS.ib.bc(KS.ps.faceCenter[i, 1], KS.ps.faceCenter[i, 2]), + face[i].n[1], + face[i].n[2], + ) KitBase.flux_boundary_maxwell!( face[i].fw, @@ -1369,7 +1343,11 @@ function evolve!( idx = ifelse(KS.ps.faceCells[i, 1] != -1, 1, 2) if KS.ps.cellType[KS.ps.faceCells[i, idx]] == 2 - bc = KitBase.local_frame(KS.ib.bc(KS.ps.faceCenter[i, 1], KS.ps.faceCenter[i, 2]), face[i].n[1], face[i].n[2]) + bc = KitBase.local_frame( + KS.ib.bc(KS.ps.faceCenter[i, 1], KS.ps.faceCenter[i, 2]), + face[i].n[1], + face[i].n[2], + ) KitBase.flux_boundary_maxwell!( face[i].fw, diff --git a/src/Solver/solver_init.jl b/src/Solver/solver_init.jl index ae44d04c..fccacca8 100644 --- a/src/Solver/solver_init.jl +++ b/src/Solver/solver_init.jl @@ -30,7 +30,7 @@ function initialize(configfilename::T) where {T<:AbstractString} else ks = SolverSet(configfilename) cftuple = init_fvm(ks) - + return (ks, cftuple..., 0.0) end end @@ -57,7 +57,11 @@ end Initialize finite volume method """ -function init_fvm(KS::T, array = :dynamic_array; structarray = false) where {T<:AbstractSolverSet} +function init_fvm( + KS::T, + array = :dynamic_array; + structarray = false, +) where {T<:AbstractSolverSet} return init_fvm(KS, KS.ps, array; structarray = structarray) end @@ -95,10 +99,7 @@ function init_fvm( end end - ctr[i] = ControlVolume1D( - funcar(w), - funcar(prim), - ) + ctr[i] = ControlVolume1D(funcar(w), funcar(prim)) end for i = 1:KS.pSpace.nx+1 @@ -116,11 +117,7 @@ function init_fvm( prim = funcprim(w, γ) f = KS.ib.ff(KS.ps.x[i]) - ctr[i] = ControlVolume1D1F( - funcar(w), - funcar(prim), - funcar(f), - ) + ctr[i] = ControlVolume1D1F(funcar(w), funcar(prim), funcar(f)) end for i = 1:KS.pSpace.nx+1 @@ -139,12 +136,7 @@ function init_fvm( prim = funcprim(w, γ) h, b = KS.ib.ff(KS.ps.x[i]) - ctr[i] = ControlVolume1D2F( - funcar(w), - funcar(prim), - funcar(h), - funcar(b), - ) + ctr[i] = ControlVolume1D2F(funcar(w), funcar(prim), funcar(h), funcar(b)) end for i = 1:KS.pSpace.nx+1 @@ -229,7 +221,7 @@ function init_fvm( array = :dynamic_array; structarray = false, ) where {T<:AbstractSolverSet,T1<:AbstractPhysicalSpace2D} - + funcar = eval(array) funcst = ifelse(structarray, StructArray, dynamic_array) funcprim = ifelse(KS.set.nSpecies == 1, conserve_prim, mixture_conserve_prim) @@ -253,16 +245,17 @@ function init_fvm( w = KS.ib.fw(KS.ps.x[i, j], KS.ps.y[i, j]) prim = funcprim(w, KS.gas.γ) - ctr[i, j] = ControlVolume2D( - funcar(w), - funcar(prim), - ) + ctr[i, j] = ControlVolume2D(funcar(w), funcar(prim)) end for j = 1:ny for i = 1:nx n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :]) - n .= ifelse(dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, n, -n) + n .= ifelse( + dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, + n, + -n, + ) a1face[i, j] = Interface2D1F( point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :]), @@ -272,20 +265,27 @@ function init_fvm( ) end n = unit_normal(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :]) - n .= ifelse(dot(n, ps.vertices[nx, j, 2, :] .- [ps.x[nx, j], ps.y[nx, j]]) >= 0, n, -n) + n .= ifelse( + dot(n, ps.vertices[nx, j, 2, :] .- [ps.x[nx, j], ps.y[nx, j]]) >= 0, + n, + -n, + ) - a1face[nx+1, j] = - Interface2D1F( - point_distance(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :]), - n[1], - n[2], - funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])), - ) + a1face[nx+1, j] = Interface2D1F( + point_distance(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :]), + n[1], + n[2], + funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])), + ) end for i = 1:nx for j = 1:ny n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :]) - n .= ifelse(dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, n, -n) + n .= ifelse( + dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, + n, + -n, + ) a2face[i, j] = Interface2D1F( point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :]), @@ -295,15 +295,18 @@ function init_fvm( ) end n = unit_normal(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :]) - n .= ifelse(dot(n, ps.vertices[i, ny, 3, :] .- [ps.x[i, ny], ps.y[i, ny]]) >= 0, n, -n) + n .= ifelse( + dot(n, ps.vertices[i, ny, 3, :] .- [ps.x[i, ny], ps.y[i, ny]]) >= 0, + n, + -n, + ) - a2face[i, ny+1] = - Interface2D1F( - point_distance(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :]), - n[1], - n[2], - funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])), - ) + a2face[i, ny+1] = Interface2D1F( + point_distance(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :]), + n[1], + n[2], + funcar(KS.ib.fw(KS.ps.x[1], KS.ps.y[1])), + ) end elseif KS.set.space[3:4] == "1f" @@ -321,17 +324,17 @@ function init_fvm( prim = funcprim(w, KS.gas.γ) h = KS.ib.ff(KS.ps.x[i, j], KS.ps.y[i, j]) - ctr[i, j] = ControlVolume2D1F( - funcar(w), - funcar(prim), - funcar(h), - ) + ctr[i, j] = ControlVolume2D1F(funcar(w), funcar(prim), funcar(h)) end for j = 1:ny for i = 1:nx n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :]) - n .= ifelse(dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, n, -n) + n .= ifelse( + dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, + n, + -n, + ) a1face[i, j] = Interface2D1F( point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :]), @@ -342,7 +345,11 @@ function init_fvm( ) end n = unit_normal(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :]) - n .= ifelse(dot(n, ps.vertices[nx, j, 2, :] .- [ps.x[nx, j], ps.y[nx, j]]) >= 0, n, -n) + n .= ifelse( + dot(n, ps.vertices[nx, j, 2, :] .- [ps.x[nx, j], ps.y[nx, j]]) >= 0, + n, + -n, + ) a1face[nx+1, j] = Interface2D1F( point_distance(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :]), @@ -355,7 +362,11 @@ function init_fvm( for i = 1:nx for j = 1:ny n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :]) - n .= ifelse(dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, n, -n) + n .= ifelse( + dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, + n, + -n, + ) a2face[i, j] = Interface2D1F( point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :]), @@ -366,7 +377,11 @@ function init_fvm( ) end n = unit_normal(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :]) - n .= ifelse(dot(n, ps.vertices[i, ny, 3, :] .- [ps.x[i, ny], ps.y[i, ny]]) >= 0, n, -n) + n .= ifelse( + dot(n, ps.vertices[i, ny, 3, :] .- [ps.x[i, ny], ps.y[i, ny]]) >= 0, + n, + -n, + ) a2face[i, ny+1] = Interface2D1F( point_distance(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :]), @@ -392,18 +407,17 @@ function init_fvm( prim = funcprim(w, KS.gas.γ) h, b = KS.ib.ff(KS.ps.x[i, j], KS.ps.y[i, j]) - ctr[i, j] = ControlVolume2D2F( - funcar(w), - funcar(prim), - funcar(h), - funcar(b), - ) + ctr[i, j] = ControlVolume2D2F(funcar(w), funcar(prim), funcar(h), funcar(b)) end for j = 1:ny for i = 1:nx n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :]) - n .= ifelse(dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, n, -n) + n .= ifelse( + dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, + n, + -n, + ) a1face[i, j] = Interface2D2F( point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 4, :]), @@ -414,7 +428,11 @@ function init_fvm( ) end n = unit_normal(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :]) - n .= ifelse(dot(n, ps.vertices[nx, j, 2, :] .- [ps.x[nx, j], ps.y[nx, j]]) >= 0, n, -n) + n .= ifelse( + dot(n, ps.vertices[nx, j, 2, :] .- [ps.x[nx, j], ps.y[nx, j]]) >= 0, + n, + -n, + ) a1face[nx+1, j] = Interface2D2F( point_distance(ps.vertices[nx, j, 2, :], ps.vertices[nx, j, 3, :]), @@ -427,7 +445,11 @@ function init_fvm( for i = 1:nx for j = 1:ny n = unit_normal(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :]) - n .= ifelse(dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, n, -n) + n .= ifelse( + dot(n, [ps.x[i, j], ps.y[i, j]] .- ps.vertices[i, j, 1, :]) >= 0, + n, + -n, + ) a2face[i, j] = Interface2D2F( point_distance(ps.vertices[i, j, 1, :], ps.vertices[i, j, 2, :]), @@ -438,7 +460,11 @@ function init_fvm( ) end n = unit_normal(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :]) - n .= ifelse(dot(n, ps.vertices[i, ny, 3, :] .- [ps.x[i, ny], ps.y[i, ny]]) >= 0, n, -n) + n .= ifelse( + dot(n, ps.vertices[i, ny, 3, :] .- [ps.x[i, ny], ps.y[i, ny]]) >= 0, + n, + -n, + ) a2face[i, ny+1] = Interface2D2F( point_distance(ps.vertices[i, ny, 3, :], ps.vertices[i, ny, 4, :]), @@ -461,7 +487,7 @@ function init_fvm( array = :dynamic_array; structarray = false, ) where {T<:AbstractSolverSet} - + funcar = eval(array) funcst = ifelse(structarray, StructArray, dynamic_array) funcprim = ifelse(KS.set.nSpecies == 1, conserve_prim, mixture_conserve_prim) @@ -717,7 +743,7 @@ function init_fvm( end return ctr |> funcst, face |> funcst - + end diff --git a/src/Solver/solver_main.jl b/src/Solver/solver_main.jl index 81d6c4fd..ea8bd6d9 100644 --- a/src/Solver/solver_main.jl +++ b/src/Solver/solver_main.jl @@ -57,7 +57,14 @@ function solve!( #dt = timestep(KS, ctr, simTime) reconstruct!(KS, ctr) - evolve!(KS, ctr, face, dt; mode = symbolize(KS.set.flux), bc = symbolize(KS.set.boundary)) + evolve!( + KS, + ctr, + face, + dt; + mode = symbolize(KS.set.flux), + bc = symbolize(KS.set.boundary), + ) update!( KS, ctr, diff --git a/src/Solver/solver_reconstruction.jl b/src/Solver/solver_reconstruction.jl index ad2d3f06..075b0e73 100644 --- a/src/Solver/solver_reconstruction.jl +++ b/src/Solver/solver_reconstruction.jl @@ -297,12 +297,7 @@ function reconstruct!( reconstruct2!(swL, ctr[1, j].w, ctr[2, j].w, 0.5 * (dx[1, j] + dx[2, j])) swR = extract_last(ctr[nx, j].sw, 1, mode = :view) - reconstruct2!( - swR, - ctr[nx-1, j].w, - ctr[nx, j].w, - 0.5 * (dx[nx-1, j] + dy[nx, j]), - ) + reconstruct2!(swR, ctr[nx-1, j].w, ctr[nx, j].w, 0.5 * (dx[nx-1, j] + dy[nx, j])) end @inbounds Threads.@threads for i = 1:nx @@ -310,12 +305,7 @@ function reconstruct!( reconstruct2!(swD, ctr[i, 1].w, ctr[i, 2].w, 0.5 * (dy[i, 1] + dy[i, 2])) swU = extract_last(ctr[i, ny].sw, 2, mode = :view) - reconstruct2!( - swU, - ctr[i, ny-1].w, - ctr[i, ny].w, - 0.5 * (dy[i, ny-1] + dy[i, ny]), - ) + reconstruct2!(swU, ctr[i, ny-1].w, ctr[i, ny].w, 0.5 * (dy[i, ny-1] + dy[i, ny])) end # inner @@ -375,12 +365,7 @@ function reconstruct!( reconstruct2!(swL, ctr[1, j].w, ctr[2, j].w, 0.5 * (dx[1, j] + dx[2, j])) swR = extract_last(ctr[nx, j].sw, 1, mode = :view) - reconstruct2!( - swR, - ctr[nx-1, j].w, - ctr[nx, j].w, - 0.5 * (dx[nx-1, j] + dx[nx, j]), - ) + reconstruct2!(swR, ctr[nx-1, j].w, ctr[nx, j].w, 0.5 * (dx[nx-1, j] + dx[nx, j])) end @inbounds Threads.@threads for i = 1:nx @@ -388,12 +373,7 @@ function reconstruct!( reconstruct2!(swD, ctr[i, 1].w, ctr[i, 2].w, 0.5 * (dy[i, 1] + dy[i, 2])) swU = extract_last(ctr[i, ny].sw, 2, mode = :view) - reconstruct2!( - swU, - ctr[i, ny-1].w, - ctr[i, ny].w, - 0.5 * (dy[i, ny-1] + dy[i, ny]), - ) + reconstruct2!(swU, ctr[i, ny-1].w, ctr[i, ny].w, 0.5 * (dy[i, ny-1] + dy[i, ny])) end # inner @@ -434,12 +414,7 @@ function reconstruct!( reconstruct2!(sfL, ctr[1, j].f, ctr[2, j].f, 0.5 * (dx[1, j] + dx[2, j])) sfR = extract_last(ctr[nx, j].sf, 1, mode = :view) - reconstruct2!( - sfR, - ctr[nx-1, j].f, - ctr[nx, j].f, - 0.5 * (dx[nx-1, j] + dx[nx, j]), - ) + reconstruct2!(sfR, ctr[nx-1, j].f, ctr[nx, j].f, 0.5 * (dx[nx-1, j] + dx[nx, j])) end @inbounds Threads.@threads for i = 1:nx @@ -447,12 +422,7 @@ function reconstruct!( reconstruct2!(sfD, ctr[i, 1].f, ctr[i, 2].f, 0.5 * (dy[i, 1] + dy[i, 2])) sfU = extract_last(ctr[i, ny].sf, 2, mode = :view) - reconstruct2!( - sfU, - ctr[i, ny-1].f, - ctr[i, ny].f, - 0.5 * (dy[i, ny-1] + dy[i, ny]), - ) + reconstruct2!(sfU, ctr[i, ny-1].f, ctr[i, ny].f, 0.5 * (dy[i, ny-1] + dy[i, ny])) end # inner @@ -512,12 +482,7 @@ function reconstruct!( reconstruct2!(swL, ctr[1, j].w, ctr[2, j].w, 0.5 * (dx[1, j] + dx[2, j])) swR = extract_last(ctr[nx, j].sw, 1, mode = :view) - reconstruct2!( - swR, - ctr[nx-1, j].w, - ctr[nx, j].w, - 0.5 * (dx[nx-1, j] + dx[nx, j]), - ) + reconstruct2!(swR, ctr[nx-1, j].w, ctr[nx, j].w, 0.5 * (dx[nx-1, j] + dx[nx, j])) end @inbounds Threads.@threads for i = 1:nx @@ -525,12 +490,7 @@ function reconstruct!( reconstruct2!(swD, ctr[i, 1].w, ctr[i, 2].w, 0.5 * (dy[i, 1] + dy[i, 2])) swU = extract_last(ctr[i, ny].sw, 2, mode = :view) - reconstruct2!( - swU, - ctr[i, ny-1].w, - ctr[i, ny].w, - 0.5 * (dy[i, ny-1] + dy[i, ny]), - ) + reconstruct2!(swU, ctr[i, ny-1].w, ctr[i, ny].w, 0.5 * (dy[i, ny-1] + dy[i, ny])) end # inner @@ -573,19 +533,9 @@ function reconstruct!( reconstruct2!(sbL, ctr[1, j].b, ctr[2, j].b, 0.5 * (dx[1, j] + dx[2, j])) shR = extract_last(ctr[nx, j].sh, 1, mode = :view) - reconstruct2!( - shR, - ctr[nx-1, j].h, - ctr[nx, j].h, - 0.5 * (dx[nx-1, j] + dx[nx, j]), - ) + reconstruct2!(shR, ctr[nx-1, j].h, ctr[nx, j].h, 0.5 * (dx[nx-1, j] + dx[nx, j])) sbR = extract_last(ctr[nx, j].sb, 1, mode = :view) - reconstruct2!( - sbR, - ctr[nx-1, j].b, - ctr[nx, j].b, - 0.5 * (dx[nx-1, j] + dx[nx, j]), - ) + reconstruct2!(sbR, ctr[nx-1, j].b, ctr[nx, j].b, 0.5 * (dx[nx-1, j] + dx[nx, j])) end @inbounds Threads.@threads for i = 1:nx @@ -595,19 +545,9 @@ function reconstruct!( reconstruct2!(sbD, ctr[i, 1].b, ctr[i, 2].b, 0.5 * (dy[i, 1] + dy[i, 2])) shU = extract_last(ctr[i, ny].sh, 2, mode = :view) - reconstruct2!( - shU, - ctr[i, ny-1].h, - ctr[i, ny].h, - 0.5 * (dy[i, ny-1] + dy[i, ny]), - ) + reconstruct2!(shU, ctr[i, ny-1].h, ctr[i, ny].h, 0.5 * (dy[i, ny-1] + dy[i, ny])) sbU = extract_last(ctr[i, ny].sb, 2, mode = :view) - reconstruct2!( - sbU, - ctr[i, ny-1].b, - ctr[i, ny].b, - 0.5 * (dy[i, ny-1] + dy[i, ny]), - ) + reconstruct2!(sbU, ctr[i, ny-1].b, ctr[i, ny].b, 0.5 * (dy[i, ny-1] + dy[i, ny])) end # inner @@ -689,12 +629,7 @@ function reconstruct!( reconstruct2!(swL, ctr[1, j].w, ctr[2, j].w, 0.5 * (dx[1, j] + dx[2, j])) swR = extract_last(ctr[nx, j].sw, 1, mode = :view) - reconstruct2!( - swR, - ctr[nx-1, j].w, - ctr[nx, j].w, - 0.5 * (dx[nx-1, j] + dx[nx, j]), - ) + reconstruct2!(swR, ctr[nx-1, j].w, ctr[nx, j].w, 0.5 * (dx[nx-1, j] + dx[nx, j])) end @inbounds Threads.@threads for i = 1:nx @@ -702,12 +637,7 @@ function reconstruct!( reconstruct2!(swD, ctr[i, 1].w, ctr[i, 2].w, 0.5 * (dy[i, 1] + dy[i, 2])) swU = extract_last(ctr[i, ny].sw, 2, mode = :view) - reconstruct2!( - swU, - ctr[i, ny-1].w, - ctr[i, ny].w, - 0.5 * (dy[i, ny-1] + dy[i, ny]), - ) + reconstruct2!(swU, ctr[i, ny-1].w, ctr[i, ny].w, 0.5 * (dy[i, ny-1] + dy[i, ny])) end # inner @@ -752,26 +682,11 @@ function reconstruct!( reconstruct2!(s2L, ctr[1, j].h2, ctr[2, j].h2, 0.5 * (dx[1, j] + dx[2, j])) s0R = extract_last(ctr[nx, j].sh0, 1, mode = :view) - reconstruct2!( - s0R, - ctr[nx-1, j].h0, - ctr[nx, j].h0, - 0.5 * (dx[nx-1, j] + dx[nx, j]), - ) + reconstruct2!(s0R, ctr[nx-1, j].h0, ctr[nx, j].h0, 0.5 * (dx[nx-1, j] + dx[nx, j])) s1R = extract_last(ctr[nx, j].sh1, 1, mode = :view) - reconstruct2!( - s1R, - ctr[nx-1, j].h1, - ctr[nx, j].h1, - 0.5 * (dx[nx-1, j] + dx[nx, j]), - ) + reconstruct2!(s1R, ctr[nx-1, j].h1, ctr[nx, j].h1, 0.5 * (dx[nx-1, j] + dx[nx, j])) s2R = extract_last(ctr[nx, j].sh2, 1, mode = :view) - reconstruct2!( - s2R, - ctr[nx-1, j].h2, - ctr[nx, j].h2, - 0.5 * (dx[nx-1, j] + dx[nx, j]), - ) + reconstruct2!(s2R, ctr[nx-1, j].h2, ctr[nx, j].h2, 0.5 * (dx[nx-1, j] + dx[nx, j])) end @inbounds Threads.@threads for i = 1:nx @@ -783,12 +698,7 @@ function reconstruct!( reconstruct2!(s2D, ctr[i, 1].h2, ctr[i, 2].h2, 0.5 * (dy[i, 1] + dy[i, 2])) s0U = extract_last(ctr[i, ny].sh0, 2, mode = :view) - reconstruct2!( - s0U, - ctr[i, ny-1].h0, - ctr[i, ny].h0, - 0.5 * (dy[i, ny-1] + dy[i, ny]), - ) + reconstruct2!(s0U, ctr[i, ny-1].h0, ctr[i, ny].h0, 0.5 * (dy[i, ny-1] + dy[i, ny])) s1U = extract_last(ctr[i, KS.pSpace.ny].sh1, 2, mode = :view) reconstruct2!( s1U, diff --git a/src/Solver/solver_set.jl b/src/Solver/solver_set.jl index 74b87369..5afe9920 100644 --- a/src/Solver/solver_set.jl +++ b/src/Solver/solver_set.jl @@ -222,7 +222,16 @@ function set_geometry(dict::T) where {T<:AbstractDict} if parse(Int, dict[:space][1]) == 1 return PSpace1D(dict[:x0], dict[:x1], dict[:nx], dict[:nxg]) elseif parse(Int, dict[:space][1]) == 2 - return PSpace2D(dict[:x0], dict[:x1], dict[:nx], dict[:y0], dict[:y1], dict[:ny], dict[:nxg], dict[:nyg]) + return PSpace2D( + dict[:x0], + dict[:x1], + dict[:nx], + dict[:y0], + dict[:y1], + dict[:ny], + dict[:nxg], + dict[:nyg], + ) else throw("No preset available for 3D simulation, please set it up manually.") end @@ -637,21 +646,21 @@ function set_ib(set, pSpace, vSpace, gas, Um = 0.15, Vm = 0.0, Tm = 1.0) "ib_" * string(set.case) end end - + ib = begin if parse(Int, set.space[3]) == 0 - fw, bc = eval(Symbol(fname))(set, pSpace, vSpace, gas) + fw, bc = config_ib(set, pSpace, vSpace, gas) IB(fw, bc) elseif parse(Int, set.space[3]) in [3, 4] && gas isa AbstractPlasma - fw, ff, fE, fB, fL, bc = eval(Symbol(fname))(set, pSpace, vSpace, gas) + fw, ff, fE, fB, fL, bc = config_ib(set, pSpace, vSpace, gas) iname = "IB" * set.space[3] * "F" eval(Symbol(iname))(fw, ff, fE, fB, fL, bc) elseif set.case == "cavity" - fw, ff, bc = eval(Symbol(fname))(set, pSpace, vSpace, gas, Um, Vm, Tm) + fw, ff, bc = config_ib(set, pSpace, vSpace, gas, Um, Vm, Tm) iname = "IB" * set.space[3] * "F" eval(Symbol(iname))(fw, ff, bc) else - fw, ff, bc = eval(Symbol(fname))(set, pSpace, vSpace, gas) + fw, ff, bc = config_ib(set, pSpace, vSpace, gas) iname = "IB" * set.space[3] * "F" eval(Symbol(iname))(fw, ff, bc) end diff --git a/src/Solver/solver_update.jl b/src/Solver/solver_update.jl index 832dce6b..e2577dde 100644 --- a/src/Solver/solver_update.jl +++ b/src/Solver/solver_update.jl @@ -390,7 +390,16 @@ function update!( end if bc isa Symbol - update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = [bc, bc], isMHD = isMHD) + update_boundary!( + KS, + ctr, + face, + dt, + residual; + coll = coll, + bc = [bc, bc], + isMHD = isMHD, + ) else update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc, isMHD = isMHD) end @@ -501,7 +510,16 @@ function update!( end if bc isa Symbol - update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = [bc, bc], isMHD = isMHD) + update_boundary!( + KS, + ctr, + face, + dt, + residual; + coll = coll, + bc = [bc, bc], + isMHD = isMHD, + ) else update_boundary!(KS, ctr, face, dt, residual; coll = coll, bc = bc, isMHD = isMHD) end @@ -563,7 +581,16 @@ function update!( end if bc isa Symbol - update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = [bc, bc, bc, bc]) + update_boundary!( + KS, + ctr, + a1face, + a2face, + dt, + residual; + coll = coll, + bc = [bc, bc, bc, bc], + ) else update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = bc) end @@ -637,7 +664,16 @@ function update!( end if bc isa Symbol - update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = [bc, bc, bc, bc]) + update_boundary!( + KS, + ctr, + a1face, + a2face, + dt, + residual; + coll = coll, + bc = [bc, bc, bc, bc], + ) else update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = bc) end @@ -718,7 +754,16 @@ function update!( end if bc isa Symbol - update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = [bc, bc, bc, bc]) + update_boundary!( + KS, + ctr, + a1face, + a2face, + dt, + residual; + coll = coll, + bc = [bc, bc, bc, bc], + ) else update_boundary!(KS, ctr, a1face, a2face, dt, residual; coll = coll, bc = bc) end diff --git a/src/Theory/theory_atom.jl b/src/Theory/theory_atom.jl index fdbaaded..fb1a9d9c 100755 --- a/src/Theory/theory_atom.jl +++ b/src/Theory/theory_atom.jl @@ -4,21 +4,22 @@ Calculate reference viscosity with variable hard sphere (VHS) model """ -ref_vhs_vis(Kn::T, alpha, omega) where {T<:Real} = +ref_vhs_vis(Kn, alpha, omega) = 5.0 * (alpha + 1.0) * (alpha + 2.0) * √π / (4.0 * alpha * (5.0 - 2.0 * omega) * (7.0 - 2.0 * omega)) * Kn """ - vhs_collision_time(prim::T, muRef, omega) where {T<:AbstractArray{<:Real,1}} + vhs_collision_time(ρ, λ, μᵣ, ω) + vhs_collision_time(prim::AbstractVector{T}, muRef, omega) where {T<:Real} Calculate collision time with variable hard sphere (VHS) model """ -vhs_collision_time(ρ::T, λ, μᵣ, ω) where {T<:Real} = μᵣ * 2.0 * λ^(1.0 - ω) / ρ +vhs_collision_time(ρ, λ, μᵣ, ω) = μᵣ * 2.0 * λ^(1.0 - ω) / ρ -vhs_collision_time(prim::T, muRef, omega) where {T<:AbstractArray{<:Real,1}} = - muRef * 2.0 * prim[end]^(1.0 - omega) / prim[1] # for rykov prim[end] should be λₜ +vhs_collision_time(prim::AbstractVector{T}, muRef, omega) where {T<:Real} = + muRef * 2.0 * prim[end]^(1.0 - omega) / prim[1] # for rykov model prim[end] should be λₜ """ @@ -33,14 +34,7 @@ vhs_collision_time(prim::T, muRef, omega) where {T<:AbstractArray{<:Real,1}} = Calculate mixture collision time from AAP model """ -function aap_hs_collision_time( - prim::T, - mi, - ni, - me, - ne, - kn, -) where {T<:AbstractArray{<:Real,2}} +function aap_hs_collision_time(prim::AbstractMatrix{T}, mi, ni, me, ne, kn) where {T<:Real} ν = similar(prim, 2) @@ -137,22 +131,22 @@ end """ - boltzmann_ode!(df, f::T, p, t) where {T<:AbstractArray{<:Real,3}} + boltzmann_ode!(df, f, p, t) RHS-ODE of Boltzmann equation """ -function boltzmann_ode!(df, f::T, p, t) where {T<:AbstractArray{<:Real,3}} +function boltzmann_ode!(df, f::AbstractArray{T,3}, p, t) where {T<:Real} Kn, M, phi, psi, phipsi = p df .= boltzmann_fft(f, Kn, M, phi, psi, phipsi) end """ - boltzmann_nuode!(df, f::T, p, t) where {T<:AbstractArray{<:Real,3}} + boltzmann_nuode!(df, f, p, t) RHS-ODE of Boltzmann equation with non-uniform velocity """ -function boltzmann_nuode!(df, f::T, p, t) where {T<:AbstractArray{<:Real,3}} +function boltzmann_nuode!(df, f::AbstractArray{T,3}, p, t) where {T<:Real} Kn, M, phi, psi, phipsi, u, v, w, vnu, u1, v1, w1, vuni = p nu = length(u) @@ -172,11 +166,11 @@ end """ - bgk_ode!(df, f::T, p, t) where {T<:AbstractArray} + bgk_ode!(df, f, p, t) RHS-ODE of BGK equation """ -function bgk_ode!(df, f::T, p, t) where {T<:AbstractArray} +function bgk_ode!(df, f::AbstractArray{T}, p, t) where {T<:Real} g, τ = p df .= (g .- f) ./ τ end diff --git a/src/Theory/theory_continuum.jl b/src/Theory/theory_continuum.jl index ac28ea94..1e09c739 100755 --- a/src/Theory/theory_continuum.jl +++ b/src/Theory/theory_continuum.jl @@ -3,19 +3,15 @@ # ============================================================ """ - prim_conserve(prim::T, γ) where {T<:AbstractArray{<:Real,1}} - + prim_conserve(prim::AbstractVector{T}, γ) where {T<:Real} prim_conserve(ρ, U, λ, γ) - prim_conserve(ρ, U, V, λ, γ) - prim_conserve(ρ, U, V, W, λ, γ) Transform primitive -> conservative variables """ -function prim_conserve(prim::T, γ) where {T<:AbstractArray{<:Real,1}} - - if eltype(prim) <: Int +function prim_conserve(prim::AbstractVector{T}, γ) where {T<:Real} + if eltype(prim) <: Integer W = similar(prim, Float64) else W = similar(prim) @@ -43,7 +39,6 @@ function prim_conserve(prim::T, γ) where {T<:AbstractArray{<:Real,1}} end return W - end prim_conserve(ρ, U, λ, γ) = prim_conserve([ρ, U, λ], γ) @@ -53,9 +48,8 @@ prim_conserve(ρ, U, V, λ, γ) = prim_conserve([ρ, U, V, λ], γ) prim_conserve(ρ, U, V, W, λ, γ) = prim_conserve([ρ, U, V, W, λ], γ) #--- Rykov ---# -function prim_conserve(prim::T, γ, Kr) where {T<:AbstractArray{<:Real,1}} - - if eltype(prim) <: Int +function prim_conserve(prim::AbstractVector{T}, γ, Kr) where {T<:Real} + if eltype(prim) <: Integer W = similar(prim, Float64, length(prim) - 1) else W = similar(prim, length(prim) - 1) @@ -77,18 +71,16 @@ function prim_conserve(prim::T, γ, Kr) where {T<:AbstractArray{<:Real,1}} end return W - end """ - mixture_prim_conserve(prim::T, γ) where {T<:AbstractArray{<:Real,2}} + mixture_prim_conserve(prim::AbstractMatrix{T}, γ) where {T<:Real} Transform multi-component primitive -> conservative variables """ -function mixture_prim_conserve(prim::T, γ) where {T<:AbstractArray{<:Real,2}} - - if eltype(prim) <: Int +function mixture_prim_conserve(prim::AbstractMatrix{T}, γ) where {T<:Real} + if eltype(prim) <: Integer w = similar(prim, Float64) else w = similar(prim) @@ -99,7 +91,6 @@ function mixture_prim_conserve(prim::T, γ) where {T<:AbstractArray{<:Real,2}} end return w - end @@ -124,9 +115,8 @@ conserve_prim(u) = [u, 0.5 * u, 1.0] conserve_prim(u, a) = [u, a, 1.0] -function conserve_prim(W::T, γ) where {T<:AbstractArray{<:Real,1}} - - if eltype(W) <: Int +function conserve_prim(W::AbstractVector{T}, γ) where {T<:Real} + if eltype(W) <: Integer prim = similar(W, Float64) else prim = similar(W) @@ -152,7 +142,6 @@ function conserve_prim(W::T, γ) where {T<:AbstractArray{<:Real,1}} end return prim - end conserve_prim(ρ, M, E, γ) = conserve_prim([ρ, M, E], γ) @@ -162,9 +151,8 @@ conserve_prim(ρ, MX, MY, E, γ) = conserve_prim([ρ, MX, MY, E], γ) conserve_prim(ρ, MX, MY, MZ, E, γ) = conserve_prim([ρ, MX, MY, MZ, E], γ) #--- Rykov ---# -function conserve_prim(w::T, K, Kr) where {T<:AbstractArray{<:Real,1}} - - if eltype(w) <: Int +function conserve_prim(w::AbstractVector{T}, K, Kr) where {T<:Real} + if eltype(w) <: Integer prim = similar(w, Float64, length(w) + 1) else prim = similar(w, length(w) + 1) @@ -188,7 +176,6 @@ function conserve_prim(w::T, K, Kr) where {T<:AbstractArray{<:Real,1}} end return prim - end @@ -197,9 +184,8 @@ end Transform multi-component conservative -> primitive variables """ -function mixture_conserve_prim(W::T, γ) where {T<:AbstractArray{<:Real,2}} - - if eltype(W) <: Int +function mixture_conserve_prim(W::AbstractMatrix{T}, γ) where {T<:Real} + if eltype(W) <: Integer prim = similar(W, Float64) else prim = similar(W) @@ -210,7 +196,6 @@ function mixture_conserve_prim(W::T, γ) where {T<:AbstractArray{<:Real,2}} end return prim - end @@ -228,14 +213,14 @@ end Calculate mixture primitive variables from AAP model """ function aap_hs_prim( - prim::X, - tau::Y, + prim::AbstractMatrix{T1}, + tau::AbstractVector{T2}, mi, ni, me, ne, kn, -) where {X<:AbstractArray{<:Real,2},Y<:AbstractArray{<:Real,1}} +) where {T1<:Real,T2<:Real} mixprim = similar(prim) @@ -465,19 +450,18 @@ burgers_flux(u) = 0.5 * u^2 """ - euler_flux(w::A, γ; frame = :cartesian::Symbol) where {A<:AbstractArray{<:Real,1}} + euler_flux(w::AbstractVector{T}, γ; frame = :cartesian::Symbol) where {T<:Real} Theoretical fluxes of Euler Equations * @return: flux tuple """ -function euler_flux(w::T, γ; frame = :cartesian::Symbol) where {T<:AbstractArray{<:Real,1}} - +function euler_flux(w::AbstractVector{T}, γ; frame = :cartesian::Symbol) where {T<:Real} prim = conserve_prim(w, γ) p = 0.5 * prim[1] / prim[end] - if eltype(w) <: Int + if eltype(w) <: Integer F = similar(w, Float64) G = similar(w, Float64) H = similar(w, Float64) @@ -526,21 +510,19 @@ function euler_flux(w::T, γ; frame = :cartesian::Symbol) where {T<:AbstractArra return F, G, H end - end """ - euler_jacobi(w::T, γ) where {T<:AbstractArray{<:Real,1}} + euler_jacobi(w::AbstractVector{T}, γ) where {T<:Real} Flux Jacobian of Euler Equations * @return: Jacobian matrix A """ -function euler_jacobi(w::T, γ) where {T<:AbstractArray{<:Real,1}} - - if eltype(w) <: Int +function euler_jacobi(w::AbstractVector{T}, γ) where {T<:Real} + if eltype(w) <: Integer A = similar(w, Float64, 3, 3) else A = similar(w, 3, 3) @@ -559,5 +541,4 @@ function euler_jacobi(w::T, γ) where {T<:AbstractArray{<:Real,1}} A[3, 3] = γ * w[2] / w[1] return A - end diff --git a/src/Theory/theory_diatomic.jl b/src/Theory/theory_diatomic.jl index 0c3b698d..70e5540b 100644 --- a/src/Theory/theory_diatomic.jl +++ b/src/Theory/theory_diatomic.jl @@ -87,10 +87,10 @@ function rykov!( ω0, ω1, ) where { - T1<:AbstractArray{<:Real,1}, - T2<:AbstractArray{<:Real,1}, - T3<:AbstractArray{<:Real,1}, - T4<:AbstractArray{<:Real,1}, + T1<:AbstractVector{<:Real}, + T2<:AbstractVector{<:Real}, + T3<:AbstractVector{<:Real}, + T4<:AbstractVector{<:Real}, } @. Ht_plus = diff --git a/src/Theory/theory_fsm.jl b/src/Theory/theory_fsm.jl index f3bca184..6ae1de07 100644 --- a/src/Theory/theory_fsm.jl +++ b/src/Theory/theory_fsm.jl @@ -216,9 +216,9 @@ end ψ::Y, phipsi::Z, ) where { - X<:AbstractArray{<:AbstractFloat,3}, - Y<:AbstractArray{<:AbstractFloat,4}, - Z<:AbstractArray{<:AbstractFloat,3}, + X<:AbstractArray{<:Real,3}, + Y<:AbstractArray{<:Real,4}, + Z<:AbstractArray{<:Real,3}, I<:Integer, } @@ -233,9 +233,9 @@ function boltzmann_fft( ψ::Y, phipsi::Z, ) where { - X<:AbstractArray{<:AbstractFloat,3}, - Y<:AbstractArray{<:AbstractFloat,4}, - Z<:AbstractArray{<:AbstractFloat,3}, + X<:AbstractArray{<:Real,3}, + Y<:AbstractArray{<:Real,4}, + Z<:AbstractArray{<:Real,3}, I<:Integer, } @@ -278,12 +278,12 @@ end ψ::TY, phipsi::TZ, ) where { - T1<:AbstractArray{<:AbstractFloat,3}, - T2<:AbstractArray{<:AbstractFloat,3}, + T1<:AbstractArray{<:Real,3}, + T2<:AbstractArray{<:Real,3}, TR<:Real, TI<:Integer, - TY<:AbstractArray{<:AbstractFloat,4}, - TZ<:AbstractArray{<:AbstractFloat,3}, + TY<:AbstractArray{<:Real,4}, + TZ<:AbstractArray{<:Real,3}, } Calculate collision operator with FFT-based fast spectral method @@ -298,12 +298,12 @@ function boltzmann_fft!( ψ::TY, phipsi::TZ, ) where { - T1<:AbstractArray{<:AbstractFloat,3}, - T2<:AbstractArray{<:AbstractFloat,3}, + T1<:AbstractArray{<:Real,3}, + T2<:AbstractArray{<:Real,3}, TR<:Real, TI<:Integer, - TY<:AbstractArray{<:AbstractFloat,4}, - TZ<:AbstractArray{<:AbstractFloat,3}, + TY<:AbstractArray{<:Real,4}, + TZ<:AbstractArray{<:Real,3}, } f_spec = f .+ 0im diff --git a/src/Theory/theory_maxwellian.jl b/src/Theory/theory_maxwellian.jl index f079593c..497ef200 100644 --- a/src/Theory/theory_maxwellian.jl +++ b/src/Theory/theory_maxwellian.jl @@ -33,21 +33,18 @@ maxwellian(u, v, w, prim::AbstractVector{T}) where {T<:Real} = """ + * Maxwellian maxwellian!(M, u, ρ, U, λ) + maxwellian!(M, u, ρ, U, V, λ) + maxwellian!(M, u, ρ, U, V, W, λ) maxwellian!(M, u, prim) - - # Rykov - maxwellian!(Ht, Bt, Rt, Hr, Br, Rr, u, prim, K, Kr) - - maxwellian!(M, u, v, ρ, U, V, λ) maxwellian!(M, u, v, prim) + maxwellian!(M, u, v, w, prim) - # Rykov + * Rykov + maxwellian!(Ht, Bt, Rt, Hr, Br, Rr, u, prim, K, Kr) maxwellian!(Ht, Bt, Rt, Hr, Br, Rr, u, v, prim, K, Kr) - maxwellian!(M, u, v, ρ, U, V, W, λ) - maxwellian!(M, u, v, w, prim) - In-place Maxwellian * @args: particle velocity quadrature points @@ -185,27 +182,16 @@ maxwellian!( """ - mixture_maxwellian(u::X, prim::Y) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:Real,2}} - - mixture_maxwellian( - u::X, - v::X, - prim::Y, - ) where {X<:AbstractArray{<:AbstractFloat,3},Y<:AbstractArray{<:Real,2}} - - mixture_maxwellian( - u::X, - v::X, - w::X, - prim::Y, - ) where {X<:AbstractArray{<:AbstractFloat,4},Y<:AbstractArray{<:Real,2}} + mixture_maxwellian(u, prim) + mixture_maxwellian(u, v, prim) + mixture_maxwellian(u, v, w, prim) Multi-component Maxwellian in discrete form """ function mixture_maxwellian( - u::X, - prim::Y, -) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:Real,2}} + u::AbstractMatrix{T1}, + prim::AbstractMatrix{T2}, +) where {T1<:Real,T2<:Real} mixM = similar(u) for j in axes(mixM, 2) @@ -219,8 +205,8 @@ end function mixture_maxwellian( u::X, v::X, - prim::Y, -) where {X<:AbstractArray{<:AbstractFloat,3},Y<:AbstractArray{<:Real,2}} + prim::AbstractMatrix{Y}, +) where {X<:AbstractArray{<:AbstractFloat,3},Y<:Real} mixM = similar(u) for k in axes(mixM, 3) @@ -235,8 +221,8 @@ function mixture_maxwellian( u::X, v::X, w::X, - prim::Y, -) where {X<:AbstractArray{<:AbstractFloat,4},Y<:AbstractArray{<:Real,2}} + prim::AbstractMatrix{Y}, +) where {X<:AbstractArray{<:AbstractFloat,4},Y<:Real} mixM = similar(u) for l in axes(mixM, 4) @@ -250,51 +236,18 @@ end """ - mixture_maxwellian!( - M::T1, - u::T2, - prim::T3, - ) where { - T1<:AbstractArray{<:AbstractFloat,2}, - T2<:AbstractArray{<:AbstractFloat,2}, - T3<:AbstractArray{<:Real,2}, - } - - mixture_maxwellian!( - M::T1, - u::T2, - v::T2, - prim::T3, - ) where { - T1<:AbstractArray{<:AbstractFloat,3}, - T2<:AbstractArray{<:AbstractFloat,3}, - T3<:AbstractArray{<:Real,2}, - } - - mixture_maxwellian!( - M::T1, - u::T2, - v::T2, - w::T2, - prim::T3, - ) where { - T1<:AbstractArray{<:AbstractFloat,4}, - T2<:AbstractArray{<:AbstractFloat,4}, - T3<:AbstractArray{<:Real,2}, - } + mixture_maxwellian!(M, u, prim) + mixture_maxwellian!(M, u, v, prim) + mixture_maxwellian!(M, u, v, w, prim) In-place multi-component Maxwellian """ function mixture_maxwellian!( - M::T1, + M::AbstractMatrix{T1}, u::T2, - prim::T3, -) where { - T1<:AbstractArray{<:AbstractFloat,2}, - T2<:AbstractArray{<:AbstractFloat,2}, - T3<:AbstractArray{<:Real,2}, -} + prim::AbstractMatrix{T3}, +) where {T1<:Real,T2<:AbstractArray{<:Real,2},T3<:Real} for j in axes(M, 2) _M = @view M[:, j] @@ -306,15 +259,11 @@ function mixture_maxwellian!( end function mixture_maxwellian!( - M::T1, + M::AbstractArray{T1,3}, u::T2, v::T2, - prim::T3, -) where { - T1<:AbstractArray{<:AbstractFloat,3}, - T2<:AbstractArray{<:AbstractFloat,3}, - T3<:AbstractArray{<:Real,2}, -} + prim::AbstractMatrix{T3}, +) where {T1<:Real,T2<:AbstractArray{<:Real,3},T3<:Real} for k in axes(M, 3) _M = @view M[:, :, k] @@ -326,16 +275,12 @@ function mixture_maxwellian!( end function mixture_maxwellian!( - M::T1, + M::AbstractArray{T1,4}, u::T2, v::T2, w::T2, - prim::T3, -) where { - T1<:AbstractArray{<:AbstractFloat,4}, - T2<:AbstractArray{<:AbstractFloat,4}, - T3<:AbstractArray{<:Real,2}, -} + prim::AbstractMatrix{T3}, +) where {T1<:Real,T2<:AbstractArray{<:Real,4},T3<:Real} for l in axes(M, 4) _M = @view M[:, :, :, l] diff --git a/src/Theory/theory_moments_mixture.jl b/src/Theory/theory_moments_mixture.jl index 09c9033e..c9dfa923 100644 --- a/src/Theory/theory_moments_mixture.jl +++ b/src/Theory/theory_moments_mixture.jl @@ -1,11 +1,10 @@ """ - mixture_gauss_moments(prim::T, inK) where {T<:AbstractArray{<:Real,2}} + mixture_gauss_moments(prim::AbstractMatrix{T}, inK) where {T<:Real} Calculate moments of Gaussian distribution in multi-component gas """ -function mixture_gauss_moments(prim::T, inK) where {T<:AbstractArray{<:Real,2}} - - if eltype(prim) <: Int +function mixture_gauss_moments(prim::AbstractMatrix{T}, inK) where {T<:Real} + if eltype(prim) <: Integer Mu = OffsetArray(similar(prim, Float64, 7, axes(prim, 2)), 0:6, axes(prim, 2)) else Mu = OffsetArray(similar(prim, 7, axes(prim, 2)), 0:6, axes(prim, 2)) @@ -62,7 +61,6 @@ function mixture_gauss_moments(prim::T, inK) where {T<:AbstractArray{<:Real,2}} return Mu, Mv, Mw, MuL, MuR end - end @@ -72,7 +70,7 @@ end Mxi::T, alpha::I, delta::I, - ) where {T<:OffsetArray{<:AbstractFloat,2},I<:Int} + ) where {T<:OffsetArray{<:Real,2},I<:Int} function mixture_moments_conserve( Mu::T, @@ -81,7 +79,7 @@ end alpha::I, beta::I, delta::I, - ) where {T<:OffsetArray{<:AbstractFloat,2},I<:Int} + ) where {T<:OffsetArray{<:Real,2},I<:Int} Calculate conservative moments of particle distribution in multi-component gas """ @@ -90,7 +88,7 @@ function mixture_moments_conserve( Mxi::T, alpha::I, delta::I, -) where {T<:OffsetArray{<:AbstractFloat,2},I<:Int} +) where {T<:OffsetArray{<:Real,2},I<:Integer} Muv = similar(Mu, 3, size(Mu, 2)) for j in axes(Muv, 2) @@ -108,7 +106,7 @@ function mixture_moments_conserve( alpha::I, beta::I, delta::I, -) where {T<:OffsetArray{<:AbstractFloat,2},I<:Int} +) where {T<:OffsetArray{<:Real,2},I<:Integer} Muv = ifelse(size(Mw, 1) == 3, similar(Mu, 4, size(Mu, 2)), similar(Mu, 5, size(Mu, 2))) for j in axes(Muv, 2) diff --git a/src/Theory/theory_moments_pure.jl b/src/Theory/theory_moments_pure.jl index 03abe113..771da2a4 100644 --- a/src/Theory/theory_moments_pure.jl +++ b/src/Theory/theory_moments_pure.jl @@ -468,116 +468,44 @@ end """ -Discrete moments of particle distribution - -* `discrete_moments(f, ω)`: direct quadrature -* `discrete_moments(f, u, ω, n)`: velocity moments + * direct quadrature + discrete_moments(f, ω) -""" -discrete_moments( - f::X, - ω::T, -) where {X<:AbstractArray{<:AbstractFloat,1},T<:AbstractArray{<:AbstractFloat,1}} = - sum(@. ω * f) + * velocity moments + discrete_moments(f, u, ω, n) -#--- 1V ---# -discrete_moments( - f::X, - u::T, - ω::T, - n, -) where {X<:AbstractArray{<:AbstractFloat,1},T<:AbstractArray{<:AbstractFloat,1}} = - sum(@. ω * u^n * f) +Discrete moments of particle distribution -#--- 2V ---# -discrete_moments( - f::X, - u::T, - ω::T, - n, -) where {X<:AbstractArray{<:AbstractFloat,2},T<:AbstractArray{<:AbstractFloat,2}} = - sum(@. ω * u^n * f) +""" +discrete_moments(f, ω) = sum(@. ω * f) -#--- 3V ---# -discrete_moments( - f::X, - u::T, - ω::T, - n, -) where {X<:AbstractArray{<:AbstractFloat,3},T<:AbstractArray{<:AbstractFloat,3}} = - sum(@. ω * u^n * f) +discrete_moments(f, u, ω, n) = sum(@. ω * u^n * f) """ + pressure(f, prim, u, ω) + pressure(h, b, prim, u, ω, K) + pressure(h, b, prim, u, v, ω, K) + Calculate pressure from particle distribution function """ -pressure( - f::X, - prim::Y, - u::Z, - ω::Z, -) where { - X<:AbstractArray{<:AbstractFloat,1}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:AbstractFloat,1}, -} = sum(@. ω * (u - prim[2])^2 * f) +pressure(f, prim, u, ω) = sum(@. ω * (u - prim[2])^2 * f) -pressure( - h::X, - b::X, - prim::Y, - u::Z, - ω::Z, - K, -) where { - X<:AbstractArray{<:AbstractFloat,1}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:AbstractFloat,1}, -} = (sum(@. ω * (u - prim[2])^2 * h) + sum(@. ω * b)) / (K + 1.0) +pressure(h, b, prim, u, ω, K) = + (sum(@. ω * (u - prim[2])^2 * h) + sum(@. ω * b)) / (K + 1.0) -pressure( - h::X, - b::X, - prim::Y, - u::Z, - v::Z, - ω::Z, - K, -) where { - X<:AbstractArray{<:AbstractFloat,2}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:AbstractFloat,2}, -} = (sum(@. ω * ((u - prim[2])^2 + (v - prim[3])^2) * h) + sum(@. ω * b)) / (K + 2.0) +pressure(h, b, prim, u, v, ω, K) = + (sum(@. ω * ((u - prim[2])^2 + (v - prim[3])^2) * h) + sum(@. ω * b)) / (K + 2.0) """ Calculate stress tensor from particle distribution function """ -stress( - f::X, - prim::Y, - u::Z, - ω::Z, -) where { - X<:AbstractArray{<:AbstractFloat,1}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:AbstractFloat,1}, -} = sum(@. ω * (u - prim[2]) * (u - prim[2]) * f) - -function stress( - f::X, - prim::Y, - u::Z, - v::Z, - ω::Z, -) where { - X<:AbstractArray{<:AbstractFloat,2}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:AbstractFloat,2}, -} +stress(f, prim, u, ω) = sum(@. ω * (u - prim[2]) * (u - prim[2]) * f) +function stress(f, prim, u, v, ω) P = similar(prim, 2, 2) P[1, 1] = sum(@. ω * (u - prim[2]) * (u - prim[2]) * f) @@ -586,24 +514,25 @@ function stress( P[1, 2] = sum(@. ω * (v - prim[3]) * (v - prim[3]) * f) return P - end """ + heat_flux(f, prim, u, ω) + heat_flux(h, b, prim, u, ω) + heat_flux(h, b, r, prim, u, ω) + heat_flux(f, prim, u, v, ω) + heat_flux(h, b, prim, u, v, ω) + heat_flux(h, b, r, prim, u, v, ω) + heat_flux(f, prim, u, v, w, ω) + Calculate heat flux from particle distribution function +Multiple dispatch doesn't consider unstructured multi-dimensional velocity space. +In that case a new method needs to be defined. + """ -heat_flux( - h::X, - prim::Y, - u::Z, - ω::Z, -) where { - X<:AbstractArray{<:AbstractFloat,1}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:AbstractFloat,1}, -} = 0.5 * sum(@. ω * (u - prim[2]) * (u - prim[2])^2 * h) # 1F1V +heat_flux(f, prim, u, ω) = 0.5 * sum(@. ω * (u - prim[2]) * (u - prim[2])^2 * f) # 1F1V #--- 2F1V ---# heat_flux( @@ -618,7 +547,7 @@ heat_flux( Z<:AbstractArray{<:AbstractFloat,1}, } = 0.5 * (sum(@. ω * (u - prim[2]) * (u - prim[2])^2 * h) + sum(@. ω * (u - prim[2]) * b)) -#--- 3F1V (Rykov) ---# +#--- 3F1V Rykov ---# function heat_flux( h::X, b::X, @@ -695,7 +624,7 @@ function heat_flux( end -#--- 3F2V (Rykov) ---# +#--- 3F2V Rykov ---# function heat_flux( h::X, b::X, diff --git a/src/Theory/theory_pdf.jl b/src/Theory/theory_pdf.jl index 7a1f31af..af853d13 100644 --- a/src/Theory/theory_pdf.jl +++ b/src/Theory/theory_pdf.jl @@ -1,17 +1,16 @@ """ pdf_slope(u, Δ) - - pdf_slope(prim::A, sw::B, inK) where {A<:AbstractArray{<:Real,1},B<:AbstractArray{<:Real,1}} + pdf_slope(prim, sw, inK) Calculate slope of particle distribution function `a = a1 + u * a2 + 0.5 * u^2 * a3` """ -pdf_slope(u, Δ::T) where {T<:Real} = Δ / (u + 1e-7) +pdf_slope(u, Δ) = Δ / (u + 1e-7) function pdf_slope( - prim::X, - sw::Y, + prim::AbstractVector{T1}, + sw::AbstractVector{T2}, inK, -) where {X<:AbstractArray{<:Real,1},Y<:AbstractArray{<:Real,1}} +) where {T1<:Real,T2<:Real} sl = similar(sw, axes(prim)) @@ -63,16 +62,16 @@ end """ - mixture_pdf_slope(prim::X, sw::Y, inK) where {X<:AbstractArray{<:Real,2},Y<:AbstractArray{<:Real,2}} + mixture_pdf_slope(prim, sw, inK) Calculate slope of multi-component particle distribution function `a = a1 + u * a2 + 0.5 * u^2 * a3` """ function mixture_pdf_slope( - prim::X, - sw::Y, + prim::AbstractMatrix{T1}, + sw::AbstractMatrix{T2}, inK, -) where {X<:AbstractArray{<:Real,2},Y<:AbstractArray{<:Real,2}} +) where {T1<:Real,T2<:Real} sl = similar(sw, axes(prim)) @@ -86,17 +85,8 @@ end """ - reduce_distribution( - f::X, - weights::Y, - dim = 1, - ) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:AbstractFloat,1}} - - reduce_distribution( - f::X, - weights::Y, - dim = 1, - ) where {X<:AbstractArray{<:AbstractFloat,3},Y<:AbstractArray{<:AbstractFloat,2}} + reduce_distribution(f, weights, dim) + reduce_distribution(f, v, w, weights, dim) Reduced distribution function @@ -105,10 +95,10 @@ Reduced distribution function """ function reduce_distribution( - f::X, - weights::Y, + f::AbstractMatrix{T1}, + weights::AbstractVector{T2}, dim = 1, -) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:AbstractFloat,1}} +) where {T1<:Real,T2<:Real} if dim == 1 h = similar(f, axes(f, 1)) @@ -129,10 +119,10 @@ function reduce_distribution( end function reduce_distribution( - f::X, - weights::Y, + f::AbstractArray{T1,3}, + weights::AbstractMatrix{T2}, dim = 1, -) where {X<:AbstractArray{<:AbstractFloat,3},Y<:AbstractArray{<:AbstractFloat,2}} +) where {T1<:Real,T2<:Real} if dim == 1 h = similar(f, axes(f, 1)) @@ -157,15 +147,15 @@ function reduce_distribution( end function reduce_distribution( - f::X, + f::AbstractArray{X,3}, v::Y, w::Y, - weights::Z, + weights::AbstractMatrix{Z}, dim = 1, ) where { - X<:AbstractArray{<:AbstractFloat,3}, - Y<:AbstractArray{<:AbstractFloat,3}, - Z<:AbstractArray{<:AbstractFloat,2}, + X<:Real, + Y<:AbstractArray{<:Real,3}, + Z<:Real, } if dim == 1 @@ -189,7 +179,6 @@ function reduce_distribution( h[k] = sum(@. weights * f[:, :, k]) b[k] = sum(@. weights * (v[:, :, k]^2 + w[:, :, k]^2) * f[:, :, k]) end - else end return h, b @@ -198,36 +187,8 @@ end """ - full_distribution( - h::X, - b::X, - u::Y, - weights::Y, - v::Z, - w::Z, - ρ, - γ = 5 / 3, - ) where { - X<:AbstractArray{<:AbstractFloat,1}, - Y<:AbstractArray{<:AbstractFloat,1}, - Z<:AbstractArray{<:AbstractFloat,3}, - } - - full_distribution( - h::X, - b::X, - u::Y, - weights::Y, - v::Z, - w::Z, - prim::A, - γ = 5 / 3, - ) where { - X<:AbstractArray{<:AbstractFloat,1}, - Y<:AbstractArray{<:AbstractFloat,1}, - Z<:AbstractArray{<:AbstractFloat,3}, - A<:AbstractArray{<:Real,1}, - } + full_distribution(h, b, u, weights, v, w, ρ, γ) + full_distribution(h, b, u, weights, v, w, prim, γ) Recover full distribution function from reduced ones @@ -287,23 +248,11 @@ full_distribution( """ - shift_pdf!( - f::T, - a, - du, - dt, - ) where {T<:AbstractArray{<:AbstractFloat,1}} - - shift_pdf!( - f::X, - a::Y, - du::Z, - dt, - ) where {X<:AbstractArray{<:AbstractFloat,2},Y<:AbstractArray{<:Real,1},Z<:AbstractArray{<:AbstractFloat,1}} + shift_pdf!(f, a, du, dt) Shift distribution function by external force """ -function shift_pdf!(f::T, a, du, dt) where {T<:AbstractArray{<:AbstractFloat,1}} +function shift_pdf!(f::AbstractVector{T}, a, du, dt) where {T<:Real} q0 = eachindex(f) |> first # for OffsetArray q1 = eachindex(f) |> last @@ -343,42 +292,25 @@ end #--- multi-component gas ---# function shift_pdf!( - f::X, - a::Y, - du::Z, + f::AbstractMatrix{X}, + a::AbstractVector{Y}, + du::AbstractVector{Z}, dt, -) where { - X<:AbstractArray{<:AbstractFloat,2}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:AbstractFloat,1}, -} +) where {X<:Real,Y<:Real,Z<:Real} + for j in axes(f, 2) _f = @view f[:, j] shift_pdf!(_f, a[j], du[j], dt) end return nothing + end """ - function chapman_enskog( - u::AbstractVector{T1}, - prim::AbstractVector{T2}, - a::AbstractVector{T3}, - A::AbstractVector{T3}, - τ::Real, - ) where {T1<:AbstractFloat,T2<:Real,T3<:Real} - - function chapman_enskog( - u::AbstractMatrix{T1}, - v::AbstractMatrix{T1}, - prim::AbstractVector{T2}, - a::AbstractVector{T3}, - b::AbstractVector{T3}, - A::AbstractVector{T3}, - τ::Real, - ) where {T1<:AbstractFloat,T2<:Real,T3<:Real} + chapman_enskog(u, prim, a, A, τ) + chapman_enskog(u, v, prim, a, b, A, τ) Recover discrete Chapman-Enskog expansion """ @@ -389,6 +321,7 @@ function chapman_enskog( A::AbstractVector{T3}, τ::Real, ) where {T1<:AbstractFloat,T2<:Real,T3<:Real} + M = maxwellian(u, prim) f = @. M * ( 1 - @@ -396,6 +329,7 @@ function chapman_enskog( ) return f + end function chapman_enskog( @@ -407,6 +341,7 @@ function chapman_enskog( A::AbstractVector{T3}, τ::Real, ) where {T1<:AbstractFloat,T2<:Real,T3<:Real} + M = maxwellian(u, v, prim) f = @. M * ( 1.0 - τ * (a[1] * u + a[2] * u^2 + a[3] * u * v + 0.5 * a[4] * u * (u^2 + v^2)) - @@ -415,4 +350,5 @@ function chapman_enskog( ) return f + end diff --git a/src/Theory/theory_plasma.jl b/src/Theory/theory_plasma.jl index 1f9aca8f..9e2ae7e5 100644 --- a/src/Theory/theory_plasma.jl +++ b/src/Theory/theory_plasma.jl @@ -1,25 +1,17 @@ """ - em_coefficients( - prim::X, - E::Y, - B::Z, - mr, - lD, - rL, - dt, - ) where {X<:AbstractArray{<:Real,2},Y<:AbstractArray{<:Real,1},Z<:AbstractArray{<:Real,1}} + em_coefficients(prim, E, B, mr, ld, rL, dt) Calculate electromagnetic coeffcients in hyperbolic Maxwell's equations """ function em_coefficients( - prim::X, - E::Y, - B::Z, + prim::AbstractMatrix{X}, + E::AbstractVector{Y}, + B::AbstractVector{Z}, mr, lD, rL, dt, -) where {X<:AbstractArray{<:Real,2},Y<:AbstractArray{<:Real,1},Z<:AbstractArray{<:Real,1}} +) where {X<:Real,Y<:Real,Z<:Real} if eltype(prim) <: Int A = zeros(9, 9) diff --git a/src/Theory/theory_shakhov.jl b/src/Theory/theory_shakhov.jl index d1b101a9..9ad18ba4 100644 --- a/src/Theory/theory_shakhov.jl +++ b/src/Theory/theory_shakhov.jl @@ -1,74 +1,9 @@ """ - shakhov( - u::X, - M::Y, - q, - prim::Z, - Pr, - ) where { - X<:AbstractArray{<:AbstractFloat,1}, - Y<:AbstractArray{<:AbstractFloat,1}, - Z<:AbstractArray{<:Real,1}, - } - - shakhov( - u::T, - H::X, - B::X, - q, - prim::Y, - Pr, - K, - ) where { - T<:AbstractArray{<:AbstractFloat,1}, - X<:AbstractArray{<:AbstractFloat,1}, - Y<:AbstractArray{<:Real,1}, - } - - shakhov( - u::T, - v::T, - M::X, - q::Y, - prim::Z, - Pr, - ) where { - T<:AbstractArray{<:AbstractFloat,2}, - X<:AbstractArray{<:AbstractFloat,2}, - Y<:AbstractArray{<:AbstractFloat,1}, - Z<:AbstractArray{<:Real,1}, - } - - shakhov( - u::T, - v::T, - H::X, - B::X, - q::Y, - prim::Z, - Pr, - K, - ) where { - T<:AbstractArray{<:AbstractFloat,2}, - X<:AbstractArray{<:AbstractFloat,2}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:Real,1}, - } - - shakhov( - u::T, - v::T, - w::T, - M::X, - q::Y, - prim::Z, - Pr, - ) where { - T<:AbstractArray{<:AbstractFloat,3}, - X<:AbstractArray{<:AbstractFloat,3}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:Real,1}, - } + shakhov(u, M, q, prim, Pr) + shakhov(u, H, B, q, prim, Pr, K) + shakhov(u, v, M, q, prim, Pr) + shakhov(u, v, H, B, q, prim, Pr, K) + shakhov(u, v, w, M, q, prim, Pr) Shakhov non-equilibrium part @@ -78,16 +13,12 @@ Shakhov non-equilibrium part """ function shakhov( - u::X, - M::Y, + u::AbstractVector{X}, + M::AbstractVector{Y}, q, - prim::Z, + prim::AbstractVector{Z}, Pr, -) where { - X<:AbstractArray{<:AbstractFloat,1}, - Y<:AbstractArray{<:AbstractFloat,1}, - Z<:AbstractArray{<:Real,1}, -} # 1F1V +) where {X<:AbstractFloat,Y<:AbstractFloat,Z<:Real} # 1F1V M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * (u - prim[2]) * @@ -101,17 +32,17 @@ end #--- 2F1V ---# function shakhov( - u::T, - H::X, - B::X, + u::AbstractVector{X}, + H::Y, + B::Y, q, - prim::Y, + prim::AbstractVector{Z}, Pr, K, ) where { - T<:AbstractArray{<:AbstractFloat,1}, - X<:AbstractArray{<:AbstractFloat,1}, - Y<:AbstractArray{<:Real,1}, + X<:AbstractFloat, + Y<:AbstractVector{<:AbstractFloat}, + Z<:Real, } H_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * @@ -133,15 +64,15 @@ end function shakhov( u::T, v::T, - M::X, - q::Y, - prim::Z, + M::AbstractMatrix{X}, + q::AbstractVector{Y}, + prim::AbstractVector{Z}, Pr, ) where { T<:AbstractArray{<:AbstractFloat,2}, - X<:AbstractArray{<:AbstractFloat,2}, - Y<:AbstractArray{<:AbstractFloat,1}, - Z<:AbstractArray{<:Real,1}, + X<:AbstractFloat, + Y<:AbstractFloat, + Z<:Real, } M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * @@ -159,15 +90,15 @@ function shakhov( v::T, H::X, B::X, - q::Y, - prim::Z, + q::AbstractVector{Y}, + prim::AbstractVector{Z}, Pr, K, ) where { T<:AbstractArray{<:AbstractFloat,2}, X<:AbstractArray{<:AbstractFloat,2}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:Real,1}, + Y<:Real, + Z<:Real, } H_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * @@ -188,15 +119,15 @@ function shakhov( u::T, v::T, w::T, - M::X, - q::Y, - prim::Z, + M::AbstractArray{X,3}, + q::AbstractVector{Y}, + prim::AbstractVector{Z}, Pr, ) where { T<:AbstractArray{<:AbstractFloat,3}, - X<:AbstractArray{<:AbstractFloat,3}, - Y<:AbstractArray{<:Real,1}, - Z<:AbstractArray{<:Real,1}, + X<:AbstractFloat, + Y<:Real, + Z<:Real, } M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * @@ -210,72 +141,11 @@ end """ -shakhov!( - S::T1, - u::T2, - M::T1, - q, - prim, - Pr, -) where {T1<:AbstractArray{<:AbstractFloat,1},T2<:AbstractArray{<:AbstractFloat,1}} - -shakhov!( - SH::T1, - SB::T1, - u::T2, - H::T1, - B::T1, - q, - prim, - Pr, - K, -) where {T1<:AbstractArray{<:AbstractFloat,1},T2<:AbstractArray{<:AbstractFloat,1}} - -shakhov!( - S::T1, - u::T2, - v::T2, - M::T1, - q::T3, - prim, - Pr, -) where { - T1<:AbstractArray{<:AbstractFloat,2}, - T2<:AbstractArray{<:AbstractFloat,2}, - T3<:AbstractArray{<:Real,1}, -} - -shakhov!( - SH::T1, - SB::T1, - u::T2, - v::T2, - H::T1, - B::T1, - q::T3, - prim, - Pr, - K, -) where { - T1<:AbstractArray{<:AbstractFloat,2}, - T2<:AbstractArray{<:AbstractFloat,2}, - T3<:AbstractArray{<:Real,1}, -} - -shakhov!( - S::T1, - u::T2, - v::T2, - w::T2, - M::T1, - q::T3, - prim, - Pr, -) where { - T1<:AbstractArray{<:AbstractFloat,3}, - T2<:AbstractArray{<:AbstractFloat,3}, - T3<:AbstractArray{<:Real,1}, -} + shakhov!(S, u, M, q, prim, Pr) + shakhov!(SH, SB, u, H, B, q, prim, Pr, K) + shakhov!(S, u, v, M, q, prim, Pr) + shakhov!(SH, SB, u, v, H, B, q, prim, Pr, K) + shakhov!(S, u, v, w, M, q, prim, Pr) In-place Shakhov non-equilibrium part @@ -286,12 +156,12 @@ In-place Shakhov non-equilibrium part """ function shakhov!( S::T1, - u::T2, + u::AbstractVector{T2}, M::T1, q, prim, Pr, -) where {T1<:AbstractArray{<:AbstractFloat,1},T2<:AbstractArray{<:AbstractFloat,1}} # 1F1V +) where {T1<:AbstractArray{<:AbstractFloat,1},T2<:AbstractFloat} # 1F1V @. S = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * (u - prim[2]) * @@ -307,14 +177,14 @@ end function shakhov!( SH::T1, SB::T1, - u::T2, + u::AbstractVector{T2}, H::T1, B::T1, q, prim, Pr, K, -) where {T1<:AbstractArray{<:AbstractFloat,1},T2<:AbstractArray{<:AbstractFloat,1}} +) where {T1<:AbstractArray{<:AbstractFloat,1},T2<:AbstractFloat} @. SH = 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * @@ -339,13 +209,13 @@ function shakhov!( u::T2, v::T2, M::T1, - q::T3, + q::AbstractVector{T3}, prim, Pr, ) where { T1<:AbstractArray{<:AbstractFloat,2}, T2<:AbstractArray{<:AbstractFloat,2}, - T3<:AbstractArray{<:Real,1}, + T3<:Real, } @. S = @@ -366,14 +236,14 @@ function shakhov!( v::T2, H::T1, B::T1, - q::T3, + q::AbstractVector{T3}, prim, Pr, K, ) where { T1<:AbstractArray{<:AbstractFloat,2}, T2<:AbstractArray{<:AbstractFloat,2}, - T3<:AbstractArray{<:Real,1}, + T3<:Real, } @. SH = @@ -398,13 +268,13 @@ function shakhov!( v::T2, w::T2, M::T1, - q::T3, + q::AbstractVector{T3}, prim, Pr, ) where { T1<:AbstractArray{<:AbstractFloat,3}, T2<:AbstractArray{<:AbstractFloat,3}, - T3<:AbstractArray{<:Real,1}, + T3<:Real, } @. S = diff --git a/src/Theory/theory_thermo.jl b/src/Theory/theory_thermo.jl index 34a4535f..61cf889c 100755 --- a/src/Theory/theory_thermo.jl +++ b/src/Theory/theory_thermo.jl @@ -3,55 +3,55 @@ # ============================================================ """ - heat_capacity_ratio(K, D::T) where {T<:Integer} + heat_capacity_ratio(K, D) + heat_capacity_ratio(K, Nr, D) Calculate heat capacity ratio """ function heat_capacity_ratio(K, D::T) where {T<:Integer} - - if D == 1 - γ = (K + 3.0) / (K + 1.0) - elseif D == 2 - γ = (K + 4.0) / (K + 2.0) - elseif D == 3 - γ = (K + 5.0) / (K + 3.0) + γ = begin + if D == 1 + (K + 3.0) / (K + 1.0) + elseif D == 2 + (K + 4.0) / (K + 2.0) + elseif D == 3 + (K + 5.0) / (K + 3.0) + end end return γ - end #--- diatomic ---# function heat_capacity_ratio(K, Nr, D::T) where {T<:Integer} - - if D == 1 - γ = (K + 3.0 + Nr) / (K + 1.0 + Nr) - elseif D == 2 - γ = (K + 4.0 + Nr) / (K + 2.0 + Nr) - elseif D == 3 - γ = (K + 5.0 + Nr) / (K + 3.0 + Nr) + γ = begin + if D == 1 + (K + 3.0 + Nr) / (K + 1.0 + Nr) + elseif D == 2 + (K + 4.0 + Nr) / (K + 2.0 + Nr) + elseif D == 3 + (K + 5.0 + Nr) / (K + 3.0 + Nr) + end end return γ - end """ - sound_speed(λ::Real, γ::Real) - sound_speed(prim::T, γ) where {T<:AbstractArray{<:Real,1}} - sound_speed(prim::T, γ) where {T<:AbstractArray{<:Real,2}} + sound_speed(λ, γ) + sound_speed(prim, γ) Calculate speed of sound """ sound_speed(λ::Real, γ::Real) = (0.5 * γ / λ)^0.5 -sound_speed(prim::T, γ) where {T<:AbstractArray{<:Real,1}} = sound_speed(prim[end], γ) +sound_speed(prim::AbstractVector{T}, γ) where {T<:Real} = sound_speed(prim[end], γ) -# mixture -function sound_speed(prim::T, γ) where {T<:AbstractArray{<:Real,2}} +#--- mixture ---# +function sound_speed(prim::AbstractMatrix{T}, γ) where {T<:Real} c = similar(prim, axes(prim, 2)) for j in eachindex(c) c[j] = sound_speed(prim[end, j], γ) diff --git a/src/Type/struct_ctr.jl b/src/Type/struct_ctr.jl index 248a1a45..1722c097 100644 --- a/src/Type/struct_ctr.jl +++ b/src/Type/struct_ctr.jl @@ -23,7 +23,7 @@ mutable struct ControlVolume1D{A,B} <: AbstractControlVolume1D sw::A end -function ControlVolume1D(W::T1, PRIM::T2) where {T1,T2} +function ControlVolume1D(W, PRIM) w = deepcopy(W) prim = deepcopy(PRIM) sw = zero(W) @@ -64,7 +64,7 @@ mutable struct ControlVolume1D1F{A,B} <: AbstractControlVolume1D sf::B end -function ControlVolume1D1F(W::T1, PRIM::T1, F::T2) where {T1,T2} +function ControlVolume1D1F(W::T, PRIM::T, F) where T w = deepcopy(W) prim = deepcopy(PRIM) sw = zero(W) @@ -130,15 +130,7 @@ function ControlVolume1D2F( sh = zero(h) sb = zero(b) - return ControlVolume1D2F{typeof(w),typeof(h)}( - w, - prim, - sw, - h, - b, - sh, - sb, - ) + return ControlVolume1D2F{typeof(w),typeof(h)}(w, prim, sw, h, b, sh, sb) end @@ -225,13 +217,7 @@ function ControlVolume1D3F( ψ = 0.0 lorenz = deepcopy(L) - return ControlVolume1D3F{ - typeof(w), - typeof(h0), - typeof(E), - typeof(ϕ), - typeof(lorenz), - }( + return ControlVolume1D3F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}( w, prim, sw, @@ -274,13 +260,7 @@ function ControlVolume1D3F( ψ = nothing lorenz = nothing - return ControlVolume1D3F{ - typeof(w), - typeof(h0), - typeof(E), - typeof(ϕ), - typeof(lorenz), - }( + return ControlVolume1D3F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}( w, prim, sw, @@ -326,13 +306,7 @@ function ControlVolume1D3F( ψ = zero(B[1, :]) lorenz = deepcopy(L) - return ControlVolume1D3F{ - typeof(w), - typeof(h0), - typeof(E), - typeof(ϕ), - typeof(lorenz), - }( + return ControlVolume1D3F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}( w, prim, sw, @@ -444,13 +418,7 @@ function ControlVolume1D4F( ψ = 0.0 lorenz = deepcopy(L) - return ControlVolume1D4F{ - typeof(w), - typeof(h0), - typeof(E), - typeof(ϕ), - typeof(lorenz), - }( + return ControlVolume1D4F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}( w, prim, sw, @@ -503,13 +471,7 @@ function ControlVolume1D4F( ψ = zero(B[1, :]) lorenz = deepcopy(L) - return ControlVolume1D4F{ - typeof(w), - typeof(h0), - typeof(E), - typeof(ϕ), - typeof(lorenz), - }( + return ControlVolume1D4F{typeof(w),typeof(h0),typeof(E),typeof(ϕ),typeof(lorenz)}( w, prim, sw, @@ -617,13 +579,7 @@ function ControlVolume2D1F(W, PRIM, F::AbstractArray) #sf = zeros(eltype(F), (axes(F)..., Base.OneTo(2))) sf = slope_array(F) - return ControlVolume2D1F{typeof(w),typeof(sw),typeof(f),typeof(sf)}( - w, - prim, - sw, - f, - sf, - ) + return ControlVolume2D1F{typeof(w),typeof(sw),typeof(f),typeof(sf)}(w, prim, sw, f, sf) end function Base.show(io::IO, ctr::ControlVolume2D1F{A,B,C,D}) where {A,B,C,D} @@ -671,6 +627,7 @@ function ControlVolume2D2F( H::AbstractArray, B::AbstractArray, ) + w = deepcopy(W) prim = deepcopy(PRIM) #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2))) @@ -692,6 +649,7 @@ function ControlVolume2D2F( sh, sb, ) + end function Base.show(io::IO, ctr::ControlVolume2D2F{A,B,C,D}) where {A,B,C,D} @@ -760,6 +718,7 @@ function ControlVolume2D3F( B0::AbstractArray, L::AbstractArray, ) + w = deepcopy(W) prim = deepcopy(PRIM) #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2))) # 2D @@ -805,6 +764,7 @@ function ControlVolume2D3F( ψ, lorenz, ) + end #--- Rykov ---# @@ -815,6 +775,7 @@ function ControlVolume2D3F( H1::AbstractArray{<:AbstractFloat,2}, H2::AbstractArray{<:AbstractFloat,2}, ) + w = deepcopy(W) prim = deepcopy(PRIM) #sw = zeros(eltype(W), (axes(W)..., Base.OneTo(2))) # 2D @@ -860,6 +821,7 @@ function ControlVolume2D3F( ψ, lorenz, ) + end function Base.show( @@ -964,7 +926,7 @@ mutable struct ControlVolumeUS1F{E,F,A,B,C,D} <: AbstractUnstructControlVolume sf::D end -function ControlVolumeUS1F(N, X, DX, W, PRIM, F::T) where {T<:AbstractArray} +function ControlVolumeUS1F(N, X, DX, W, PRIM, F::AbstractArray{T}) where T n = deepcopy(N) x = deepcopy(X) dx = deepcopy(DX) @@ -1042,6 +1004,7 @@ function ControlVolumeUS2F( H::T2, B::T2, ) where {T1<:AbstractArray,T2<:AbstractArray} + n = deepcopy(N) x = deepcopy(X) dx = deepcopy(DX) @@ -1082,4 +1045,5 @@ function ControlVolumeUS2F( sh, sb, ) + end diff --git a/src/Type/struct_general.jl b/src/Type/struct_general.jl index 11b5652e..2b6f81cb 100644 --- a/src/Type/struct_general.jl +++ b/src/Type/struct_general.jl @@ -55,13 +55,7 @@ function Setup( end end - return Setup{ - typeof(matter), - typeof(ns), - typeof(boundary), - typeof(cfl), - typeof(time), - }( + return Setup{typeof(matter),typeof(ns),typeof(boundary),typeof(cfl),typeof(time)}( matter, case, space, diff --git a/src/Type/struct_ib.jl b/src/Type/struct_ib.jl index ee878d11..8333315a 100644 --- a/src/Type/struct_ib.jl +++ b/src/Type/struct_ib.jl @@ -19,7 +19,7 @@ end function IB(fw, gas::Scalar) γ = gas.a - bc = function(args...) + bc = function (args...) w = fw(args...) return ifelse(γ == 0, conserve_prim(w), conserve_prim(w, γ)) end @@ -28,7 +28,7 @@ function IB(fw, gas::Scalar) end function IB(fw, gas::AbstractGas) - bc = function(args...) + bc = function (args...) w = fw(args...) prim = begin if ndims(w) == 1 @@ -62,7 +62,7 @@ mutable struct IB1F{T} <: AbstractCondition end function IB1F(fw, vs::AbstractVelocitySpace, gas::AbstractProperty) - bc = function(args...) + bc = function (args...) w = fw(args...) prim = begin if ndims(w) == 1 @@ -75,7 +75,7 @@ function IB1F(fw, vs::AbstractVelocitySpace, gas::AbstractProperty) return prim end - ff = function(args...) + ff = function (args...) prim = bc(args...) M = begin @@ -124,7 +124,7 @@ mutable struct IB2F{T} <: AbstractCondition end function IB2F(fw, vs::AbstractVelocitySpace, gas::AbstractProperty) - bc = function(args...) + bc = function (args...) w = fw(args...) prim = begin if ndims(w) == 1 @@ -137,7 +137,7 @@ function IB2F(fw, vs::AbstractVelocitySpace, gas::AbstractProperty) return prim end - ff = function(args...) + ff = function (args...) prim = bc(args...) if !isdefined(vs, :v) diff --git a/test/test_config.jl b/test/test_config.jl index 6e0e902e..deba5121 100644 --- a/test/test_config.jl +++ b/test/test_config.jl @@ -3,4 +3,4 @@ u = collect(-5:0.1:5) KitBase.ib_briowu(5 / 3, 1.0, 0.5, hcat(u, u)) KitBase.ib_briowu(5 / 3, 1.0, 0.5, ones(8, 8, 2), ones(8, 8, 2)) -ib_rh(1.2, 5/3) +ib_rh(1.2, 5 / 3) diff --git a/test/test_io.jl b/test/test_io.jl index b54d85de..a0a7fdb2 100644 --- a/test/test_io.jl +++ b/test/test_io.jl @@ -10,4 +10,4 @@ KitBase.write_jld(ks, ctr) using KitBase.Plots plot_line(ks, ctr) -#plot(ks, ctr, legend=:none) \ No newline at end of file +#plot(ks, ctr, legend=:none) diff --git a/test/test_particle.jl b/test/test_particle.jl index 80fdacbc..952f0673 100644 --- a/test/test_particle.jl +++ b/test/test_particle.jl @@ -100,7 +100,7 @@ begin wL = KitBase.prim_conserve(primL, γ) wR = KitBase.prim_conserve(primR, γ) - fw = function(x) + fw = function (x) if x <= (pSpace.x0 + pSpace.x1) / 2 wL else diff --git a/test/test_reconstruction.jl b/test/test_reconstruction.jl index e6889b75..e1aa4baa 100644 --- a/test/test_reconstruction.jl +++ b/test/test_reconstruction.jl @@ -70,7 +70,7 @@ let deg = 2, u = rand(deg + 1) modal_filter!(u, 10; filter = :houli) end -let deg = 2, u = rand(deg+1, deg+1) +let deg = 2, u = rand(deg + 1, deg + 1) modal_filter!(u, 1e-6, 1e-6; filter = :l2) modal_filter!(u, 1e-6, 1e-6; filter = :l2opt) end diff --git a/test/test_solver_scalar.jl b/test/test_solver_scalar.jl index 51c57795..d181f2b5 100644 --- a/test/test_solver_scalar.jl +++ b/test/test_solver_scalar.jl @@ -17,7 +17,7 @@ vSpace = nothing property = Scalar(1.0, 1e-6) w0 = 1.0 prim0 = conserve_prim(w0, property.a) -ib = IB(x->sin(2π * x), property) +ib = IB(x -> sin(2π * x), property) ks = SolverSet(set, pSpace, vSpace, property, ib) ctr, face = init_fvm(ks, ks.ps) diff --git a/test/test_struct.jl b/test/test_struct.jl index 0768e460..1ce3df6a 100644 --- a/test/test_struct.jl +++ b/test/test_struct.jl @@ -151,24 +151,9 @@ KitBase.ControlVolume1D4F( KitBase.ControlVolume2D(w, prim) |> show KitBase.ControlVolume2D1F(w, prim, h) |> show KitBase.ControlVolume2D2F(w, prim, h, b) |> show -KitBase.ControlVolume2D3F( - w, - prim, - h, - b, - b, - zeros(3), - zeros(3), - zeros(3, 2), -) |> show +KitBase.ControlVolume2D3F(w, prim, h, b, b, zeros(3), zeros(3), zeros(3, 2)) |> show # Rykov -KitBase.ControlVolume2D3F( - rand(4), - rand(4), - rand(8, 8), - rand(8, 8), - rand(8, 8), -) |> show +KitBase.ControlVolume2D3F(rand(4), rand(4), rand(8, 8), rand(8, 8), rand(8, 8)) |> show # unstructured KitBase.ControlVolumeUS([1, 0], x, dx, w, prim) KitBase.ControlVolumeUS1F([1, 0], x, dx, w, prim, h) diff --git a/test/test_unstructure.jl b/test/test_unstructure.jl index 144087fe..1e9bd5ae 100644 --- a/test/test_unstructure.jl +++ b/test/test_unstructure.jl @@ -11,7 +11,7 @@ begin hL = KitBase.maxwellian(vs.u, vs.v, primL) bL = @. hL * gas.K / 2 / primL[end] - fw = function(x, y) + fw = function (x, y) return wL end ib = KitBase.IB2F(fw, vs, gas) @@ -132,7 +132,7 @@ begin primL = [1.0, KitBase.sound_speed(1.0, gas.γ) * gas.Ma, 0.0, 1.0] wL = KitBase.prim_conserve(primL, gas.γ) hL = KitBase.maxwellian(vs.u, vs.v, primL) - fw = function(x, y) + fw = function (x, y) return wL end ib = KitBase.IB1F(fw, vs, gas)