From 7ba3a287e814bc4555be943276953cc71c10b691 Mon Sep 17 00:00:00 2001 From: vavrines Date: Tue, 2 Aug 2022 23:05:24 +0200 Subject: [PATCH 1/5] Add kwargs in VSpace constructors --- src/Phase/velocity.jl | 182 +++++++++++++++++++-------------------- src/Solver/solver_set.jl | 68 +++++++-------- test/test_particle.jl | 2 +- test/test_quadrature.jl | 20 ++--- test/test_theory.jl | 2 +- 5 files changed, 137 insertions(+), 137 deletions(-) diff --git a/src/Phase/velocity.jl b/src/Phase/velocity.jl index fd191483..ee7146bb 100644 --- a/src/Phase/velocity.jl +++ b/src/Phase/velocity.jl @@ -23,36 +23,36 @@ end function VSpace1D( U0, U1, - NU::TI, - TYPE = "rectangle", - NG = zero(NU)::TI, - PRECISION = Float64, + NU::TI; + type = "rectangle", + ng = zero(NU)::TI, + precision = Float64, ) where {TI<:Integer} δ = (U1 - U0) / NU u = begin - if NG > 0 - OffsetArray{PRECISION}(undef, 1-NG:NU+NG) + if ng > 0 + OffsetArray{precision}(undef, 1-ng:NU+ng) else - Array{PRECISION}(undef, NU) + Array{precision}(undef, NU) end end du = similar(u) weights = similar(u) - if TYPE == "rectangle" # rectangular + if type == "rectangle" # rectangular for i in eachindex(u) u[i] = U0 + (i - 0.5) * δ du[i] = δ weights[i] = δ end - elseif TYPE == "newton" # newton-cotes + elseif type == "newton" # newton-cotes for i in eachindex(u) u[i] = U0 + (i - 0.5) * δ du[i] = δ - weights[i] = newton_cotes(i + NG, NU + NG * 2) * δ + weights[i] = newton_cotes(i + ng, NU + ng * 2) * δ end - elseif TYPE == "algebra" # algebraic + elseif type == "algebra" # algebraic _nu = NU + 1 _u = [U1 / (_nu - 1)^3 * (-_nu + 1 + 2 * (i - 1))^3 for i = 1:_nu] u .= (_u[1:end-1] .+ _u[2:end]) ./ 2 @@ -62,7 +62,7 @@ function VSpace1D( throw("No velocity quadrature available") end - return VSpace1D{PRECISION,TI,typeof(u),typeof(weights)}(U0, U1, NU, u, du, weights) + return VSpace1D{precision,TI,typeof(u),typeof(weights)}(U0, U1, NU, u, du, weights) end @@ -78,36 +78,36 @@ function MVSpace1D( Ui1, Ue0, Ue1, - NU::TI, - TYPE = "rectangle", - NG = zero(NU)::TI, - PRECISION = Float64, + NU::TI; + type = "rectangle", + ng = zero(NU)::TI, + precision = Float64, ) where {TI<:Integer} - u0 = PRECISION.([Ui0, Ue0]) - u1 = PRECISION.([Ui1, Ue1]) + u0 = precision.([Ui0, Ue0]) + u1 = precision.([Ui1, Ue1]) δ = (u1 .- u0) ./ NU u = begin - if NG > 0 - OffsetArray{PRECISION}(undef, 1-NG:NU+NG, 1:2) + if ng > 0 + OffsetArray{precision}(undef, 1-ng:NU+ng, 1:2) else - Array{PRECISION}(undef, NU, 2) + Array{precision}(undef, NU, 2) end end du = similar(u) weights = similar(u) - if TYPE == "rectangle" # rectangular + if type == "rectangle" # rectangular for j in axes(u, 2), i in axes(u, 1) u[i, j] = u0[j] + (i - 0.5) * δ[j] du[i, j] = δ[j] weights[i, j] = δ[j] end - elseif TYPE == "newton" # newton-cotes + elseif type == "newton" # newton-cotes for j in axes(u, 2), i in axes(u, 1) u[i, j] = u0[j] + (i - 0.5) * δ[j] du[i, j] = δ[j] - weights[i, j] = newton_cotes(i + NG, NU + NG * 2) * δ[j] + weights[i, j] = newton_cotes(i + ng, NU + ng * 2) * δ[j] end else throw("No velocity quadrature available") @@ -150,20 +150,20 @@ function VSpace2D( NU::TI, V0, V1, - NV::TI, - TYPE = "rectangle", - NGU = zero(NU)::TI, - NGV = zero(NV)::TI, - PRECISION = Float64, + NV::TI; + type = "rectangle", + ngu = zero(NU)::TI, + ngv = zero(NV)::TI, + precision = Float64, ) where {TI<:Integer} δu = (U1 - U0) / NU δv = (V1 - V0) / NV u = begin - if NGU > 0 || NGV > 0 - OffsetArray{PRECISION}(undef, 1-NGU:NU+NGU, 1-NGV:NV+NGV) + if ngu > 0 || ngv > 0 + OffsetArray{precision}(undef, 1-ngu:NU+ngu, 1-ngv:NV+ngv) else - Array{PRECISION}(undef, NU, NV) + Array{precision}(undef, NU, NV) end end v = similar(u) @@ -171,7 +171,7 @@ function VSpace2D( dv = similar(u) weights = similar(u) - if TYPE == "rectangle" # rectangular + if type == "rectangle" # rectangular for j in axes(u, 2) for i in axes(u, 1) u[i, j] = U0 + (i - 0.5) * δu @@ -181,7 +181,7 @@ function VSpace2D( weights[i, j] = δu * δv end end - elseif TYPE == "newton" # newton-cotes + elseif type == "newton" # newton-cotes for j in axes(u, 2) for i in axes(u, 1) u[i, j] = U0 + (i - 0.5) * δu @@ -189,13 +189,13 @@ function VSpace2D( du[i, j] = δu dv[i, j] = δv weights[i, j] = - newton_cotes(i + NGU, NU + NGU * 2) * + newton_cotes(i + ngu, NU + ngu * 2) * δu * - newton_cotes(j + NGV, NV + NGV * 2) * + newton_cotes(j + ngv, NV + ngv * 2) * δv end end - elseif TYPE == "algebra" + elseif type == "algebra" _nu = NU + 1 _nv = NV + 1 _u = [U1 / (_nu - 1)^3 * (-_nu + 1 + 2 * (i - 1))^3 for i = 1:_nu] @@ -223,7 +223,7 @@ function VSpace2D( #weights[i, j] = wx[i] * wy[j] end end - elseif TYPE == "maxwell" + elseif type == "maxwell" _u, _uw = maxwell_quadrature(NU, U1 / 5.76) _v, _vw = maxwell_quadrature(NV, V1 / 5.76) u1, v1 = meshgrid(_u, _v) @@ -260,7 +260,7 @@ function VSpace2D( throw("No velocity quadrature available") end - return VSpace2D{PRECISION,TI,typeof(u)}(U0, U1, NU, V0, V1, NV, u, v, du, dv, weights) + return VSpace2D{precision,TI,typeof(u)}(U0, U1, NU, V0, V1, NV, u, v, du, dv, weights) end @@ -281,24 +281,24 @@ function MVSpace2D( Vi1, Ve0, Ve1, - NV::TI, - TYPE = "rectangle", - NGU = zero(NU)::TI, - NGV = zero(NV)::TI, - PRECISION = Float64, + NV::TI; + type = "rectangle", + ngu = zero(NU)::TI, + ngv = zero(NV)::TI, + precision = Float64, ) where {TI<:Integer} - u0 = PRECISION.([Ui0, Ue0]) - u1 = PRECISION.([Ui1, Ue1]) + u0 = precision.([Ui0, Ue0]) + u1 = precision.([Ui1, Ue1]) δu = (u1 .- u0) ./ NU - v0 = PRECISION.([Vi0, Ve0]) - v1 = PRECISION.([Vi1, Ve1]) + v0 = precision.([Vi0, Ve0]) + v1 = precision.([Vi1, Ve1]) δv = (v1 .- v0) ./ NV u = begin - if NGU > 0 || NGV > 0 - OffsetArray{PRECISION}(undef, 1-NGU:NU+NGU, 1-NGV:NV+NGV, 1:2) + if ngu > 0 || ngv > 0 + OffsetArray{precision}(undef, 1-ngu:NU+ngu, 1-ngv:NV+ngv, 1:2) else - Array{PRECISION}(undef, NU, NV, 2) + Array{precision}(undef, NU, NV, 2) end end v = similar(u) @@ -306,7 +306,7 @@ function MVSpace2D( dv = similar(u) weights = similar(u) - if TYPE == "rectangle" # rectangular + if type == "rectangle" # rectangular for k in axes(u, 3), j in axes(u, 2), i in axes(u, 1) u[i, j, k] = u0[k] + (i - 0.5) * δu[k] v[i, j, k] = v0[k] + (j - 0.5) * δv[k] @@ -314,16 +314,16 @@ function MVSpace2D( dv[i, j, k] = δv[k] weights[i, j, k] = δu[k] * δv[k] end - elseif TYPE == "newton" # newton-cotes + elseif type == "newton" # newton-cotes for k in axes(u, 3), j in axes(u, 2), i in axes(u, 1) u[i, j, k] = u0[k] + (i - 0.5) * δu[k] v[i, j, k] = v0[k] + (j - 0.5) * δv[k] du[i, j, k] = δu[k] dv[i, j, k] = δv[k] weights[i, j, k] = - newton_cotes(i + NGU, NU + NGU * 2) * + newton_cotes(i + ngu, NU + ngu * 2) * δu[k] * - newton_cotes(j + NGV, NV + NGV * 2) * + newton_cotes(j + ngv, NV + ngv * 2) * δv[k] end else @@ -376,22 +376,22 @@ function VSpace3D( NV::TI, W0, W1, - NW::TI, - TYPE = "rectangle", - NGU = zero(NU)::TI, - NGV = zero(NV)::TI, - NGW = zero(NW)::TI, - PRECISION = Float64, + NW::TI; + type = "rectangle", + ngu = zero(NU)::TI, + ngv = zero(NV)::TI, + ngw = zero(NW)::TI, + precision = Float64, ) where {TI<:Integer} δu = (U1 - U0) / NU δv = (V1 - V0) / NV δw = (W1 - W0) / NW u = begin - if NGU > 0 || NGV > 0 || NGW > 0 - OffsetArray{PRECISION}(undef, 1-NGU:NU+NGU, 1-NGV:NV+NGV, 1-NGW:NW+NGW) + if ngu > 0 || ngv > 0 || ngw > 0 + OffsetArray{precision}(undef, 1-ngu:NU+ngu, 1-ngv:NV+ngv, 1-ngw:NW+ngw) else - Array{PRECISION}(undef, NU, NV, NW) + Array{precision}(undef, NU, NV, NW) end end v = similar(u) @@ -401,7 +401,7 @@ function VSpace3D( dw = similar(u) weights = similar(u) - if TYPE == "rectangle" # rectangular + if type == "rectangle" # rectangular for k in axes(u, 3), j in axes(u, 2), i in axes(u, 1) u[i, j, k] = U0 + (i - 0.5) * δu v[i, j, k] = V0 + (j - 0.5) * δv @@ -411,7 +411,7 @@ function VSpace3D( dw[i, j, k] = δw weights[i, j, k] = δu * δv * δw end - elseif TYPE == "newton" # newton-cotes + elseif type == "newton" # newton-cotes for k in axes(u, 3), j in axes(u, 2), i in axes(u, 1) u[i, j, k] = U0 + (i - 0.5) * δu v[i, j, k] = V0 + (j - 0.5) * δv @@ -420,14 +420,14 @@ function VSpace3D( dv[i, j, k] = δv dw[i, j, k] = δw weights[i, j, k] = - newton_cotes(i + NGU, NU + NGU * 2) * + newton_cotes(i + ngu, NU + ngu * 2) * δu * - newton_cotes(j + NGV, NV + NGV * 2) * + newton_cotes(j + ngv, NV + ngv * 2) * δv * - newton_cotes(k + NGW, NW + NGW * 2) * + newton_cotes(k + ngw, NW + ngw * 2) * δw end - elseif TYPE == "algebra" + elseif type == "algebra" _nu = NU + 1 _nv = NV + 1 _nw = NW + 1 @@ -457,7 +457,7 @@ function VSpace3D( throw("No velocity quadrature available") end - return VSpace3D{PRECISION,TI,typeof(u)}( + return VSpace3D{precision,TI,typeof(u)}( U0, U1, NU, @@ -501,29 +501,29 @@ function MVSpace3D( Wi1, We0, We1, - NW::TI, - TYPE = "rectangle", - NGU = zero(NU)::TI, - NGV = zero(NV)::TI, - NGW = zero(NW)::TI, - PRECISION = Float64, + NW::TI; + type = "rectangle", + ngu = zero(NU)::TI, + ngv = zero(NV)::TI, + ngw = zero(NW)::TI, + precision = Float64, ) where {TI<:Integer} - u0 = PRECISION.([Ui0, Ue0]) - u1 = PRECISION.([Ui1, Ue1]) + u0 = precision.([Ui0, Ue0]) + u1 = precision.([Ui1, Ue1]) δu = (u1 .- u0) ./ NU - v0 = PRECISION.([Vi0, Ve0]) - v1 = PRECISION.([Vi1, Ve1]) + v0 = precision.([Vi0, Ve0]) + v1 = precision.([Vi1, Ve1]) δv = (v1 .- v0) ./ NV - w0 = PRECISION.([Wi0, We0]) - w1 = PRECISION.([Wi1, We1]) + w0 = precision.([Wi0, We0]) + w1 = precision.([Wi1, We1]) δw = (w1 .- w0) ./ NW u = begin - if NGU > 0 || NGV > 0 || NGW > 0 - OffsetArray{PRECISION}(undef, 1-NGU:NU+NGU, 1-NGV:NV+NGV, 1-NGW:NW+NGW, 1:2) + if ngu > 0 || ngv > 0 || ngw > 0 + OffsetArray{precision}(undef, 1-ngu:NU+ngu, 1-ngv:NV+ngv, 1-ngw:NW+ngw, 1:2) else - Array{PRECISION}(undef, NU, NV, NW, 2) + Array{precision}(undef, NU, NV, NW, 2) end end v = similar(u) @@ -533,7 +533,7 @@ function MVSpace3D( dw = similar(u) weights = similar(u) - if TYPE == "rectangle" # rectangular + if type == "rectangle" # rectangular for l in axes(u, 4), k in axes(u, 3), j in axes(u, 2), i in axes(u, 1) u[i, j, k, l] = u0[l] + (i - 0.5) * δu[l] v[i, j, k, l] = v0[l] + (j - 0.5) * δv[l] @@ -543,7 +543,7 @@ function MVSpace3D( dw[i, j, k, l] = δw[l] weights[i, j, k, l] = δu[l] * δv[l] * δw[l] end - elseif TYPE == "newton" # newton-cotes + elseif type == "newton" # newton-cotes for l in axes(u, 4), k in axes(u, 3), j in axes(u, 2), i in axes(u, 1) u[i, j, k, l] = u0[l] + (i - 0.5) * δu[l] v[i, j, k, l] = v0[l] + (j - 0.5) * δv[l] @@ -552,11 +552,11 @@ function MVSpace3D( dv[i, j, k, l] = δv[l] dw[i, j, k, l] = δw[l] weights[i, j, k, l] = - newton_cotes(i + NGU, NU + NGU * 2) * + newton_cotes(i + ngu, NU + ngu * 2) * δu[l] * - newton_cotes(j + NGV, NV + NGV * 2) * + newton_cotes(j + ngv, NV + ngv * 2) * δv[l] * - newton_cotes(k + NGW, NW + NGW * 2) * + newton_cotes(k + ngw, NW + ngw * 2) * δw[l] end else diff --git a/src/Solver/solver_set.jl b/src/Solver/solver_set.jl index ecf48632..a4968f1f 100644 --- a/src/Solver/solver_set.jl +++ b/src/Solver/solver_set.jl @@ -244,17 +244,17 @@ function set_velocity(dict::T) where {T<:AbstractDict} vSpace = nothing elseif Dv == 1 if nSpecies == 1 - vSpace = VSpace1D(umin, umax, nu, vMeshType, nug) + vSpace = VSpace1D(umin, umax, nu; type = vMeshType, ng = nug) elseif nSpecies == 2 ue0 = umin * sqrt(mi / me) ue1 = umax * sqrt(mi / me) - vSpace = MVSpace1D(umin, umax, ue0, ue1, nu, vMeshType, nug) + vSpace = MVSpace1D(umin, umax, ue0, ue1, nu; type = vMeshType, ng = nug) else throw("The velocity space only supports up to two species.") end elseif Dv == 2 if nSpecies == 1 - vSpace = VSpace2D(umin, umax, nu, vmin, vmax, nv, vMeshType, nug, nvg) + vSpace = VSpace2D(umin, umax, nu, vmin, vmax, nv; type = vMeshType, ngu = nug, ngv = nvg) elseif nSpecies == 2 ue0 = umin * sqrt(mi / me) ue1 = umax * sqrt(mi / me) @@ -270,10 +270,10 @@ function set_velocity(dict::T) where {T<:AbstractDict} vmax, ve0, ve1, - nv, - vMeshType, - nug, - nvg, + nv; + type = vMeshType, + ngu = nug, + ngv = nvg, ) else throw("The velocity space only supports up to two species.") @@ -289,11 +289,11 @@ function set_velocity(dict::T) where {T<:AbstractDict} nv, wmin, wmax, - nw, - vMeshType, - nug, - nvg, - nwg, + nw; + type = vMeshType, + ngu = nug, + ngv = nvg, + ngw = nwg, ) elseif nSpecies == 2 ue0 = umin * sqrt(mi / me) @@ -317,11 +317,11 @@ function set_velocity(dict::T) where {T<:AbstractDict} wmax, we0, we1, - nw, - vMeshType, - nug, - nvg, - nwg, + nw; + type = vMeshType, + ngu = nug, + ngv = nvg, + ngw = nwg, ) else throw("The velocity space only supports up to two species.") @@ -359,17 +359,17 @@ function set_velocity(; vSpace = nothing elseif Dv == 1 if nSpecies == 1 - vSpace = VSpace1D(umin, umax, nu, vMeshType, nug) + vSpace = VSpace1D(umin, umax, nu; type = vMeshType, ng = nug) elseif nSpecies == 2 ue0 = umin * sqrt(mi / me) ue1 = umax * sqrt(mi / me) - vSpace = MVSpace1D(umin, umax, ue0, ue1, nu, vMeshType, nug) + vSpace = MVSpace1D(umin, umax, ue0, ue1, nu; type = vMeshType, ng = nug) else throw("The velocity space only supports up to two species.") end elseif Dv == 2 if nSpecies == 1 - vSpace = VSpace2D(umin, umax, nu, vmin, vmax, nv, vMeshType, nug, nvg) + vSpace = VSpace2D(umin, umax, nu, vmin, vmax, nv; type = vMeshType, ngu = nug, ngv = nvg) elseif nSpecies == 2 ue0 = umin * sqrt(mi / me) ue1 = umax * sqrt(mi / me) @@ -385,10 +385,10 @@ function set_velocity(; vmax, ve0, ve1, - nv, - vMeshType, - nug, - nvg, + nv; + type = vMeshType, + ngu = nug, + ngv = nvg, ) else throw("The velocity space only supports up to two species.") @@ -404,11 +404,11 @@ function set_velocity(; nv, wmin, wmax, - nw, - vMeshType, - nug, - nvg, - nwg, + nw; + type = vMeshType, + ngu = nug, + ngv = nvg, + ngw = nwg, ) elseif nSpecies == 2 ue0 = umin * sqrt(mi / me) @@ -432,11 +432,11 @@ function set_velocity(; wmax, we0, we1, - nw, - vMeshType, - nug, - nvg, - nwg, + nw; + type = vMeshType, + ngu = nug, + ngv = nvg, + ngw = nwg, ) else throw("The velocity space only supports up to two species.") diff --git a/test/test_particle.jl b/test/test_particle.jl index 0b0a43a2..e0d26c43 100644 --- a/test/test_particle.jl +++ b/test/test_particle.jl @@ -90,7 +90,7 @@ begin maxTime, ) pSpace = KitBase.PSpace1D(x0, x1, nx, nxg) - vSpace = KitBase.VSpace1D(umin, umax, nu, vMeshType, nug) + vSpace = KitBase.VSpace1D(umin, umax, nu; type = vMeshType, ng = nug) μᵣ = KitBase.ref_vhs_vis(knudsen, alphaRef, omegaRef) gas = KitBase.Gas( knudsen, diff --git a/test/test_quadrature.jl b/test/test_quadrature.jl index 6f3d8793..ad5cf509 100644 --- a/test/test_quadrature.jl +++ b/test/test_quadrature.jl @@ -9,14 +9,14 @@ KitBase.MVSpace1D() |> show KitBase.MVSpace2D() |> show KitBase.MVSpace3D() |> show -KitBase.VSpace1D(-5, 5, 16, "newton") -KitBase.VSpace1D(-5, 5, 16, "algebra") -KitBase.VSpace2D(-5, 5, 16, -5, 5, 16, "newton") -KitBase.VSpace2D(-5, 5, 16, -5, 5, 16, "algebra") -KitBase.VSpace2D(-5, 5, 16, -5, 5, 16, "maxwell") -KitBase.VSpace3D(-5, 5, 16, -5, 5, 16, -5, 5, 16, "newton") -KitBase.VSpace3D(-5, 5, 16, -5, 5, 16, -5, 5, 16, "algebra") -KitBase.MVSpace1D(-5, 5, -5, 5, 16, "newton") -KitBase.MVSpace2D(-5, 5, -5, 5, 16, -5, 5, -5, 5, 16, "newton") -KitBase.MVSpace3D(-5, 5, -5, 5, 8, -5, 5, -5, 5, 8, -5, 5, -5, 5, 8, "newton") +KitBase.VSpace1D(-5, 5, 16, type = "newton") +KitBase.VSpace1D(-5, 5, 16, type = "algebra") +KitBase.VSpace2D(-5, 5, 16, -5, 5, 16, type = "newton") +KitBase.VSpace2D(-5, 5, 16, -5, 5, 16, type = "algebra") +KitBase.VSpace2D(-5, 5, 16, -5, 5, 16, type = "maxwell") +KitBase.VSpace3D(-5, 5, 16, -5, 5, 16, -5, 5, 16, type = "newton") +KitBase.VSpace3D(-5, 5, 16, -5, 5, 16, -5, 5, 16, type = "algebra") +KitBase.MVSpace1D(-5, 5, -5, 5, 16, type = "newton") +KitBase.MVSpace2D(-5, 5, -5, 5, 16, -5, 5, -5, 5, 16, type = "newton") +KitBase.MVSpace3D(-5, 5, -5, 5, 8, -5, 5, -5, 5, 8, -5, 5, -5, 5, 8, type = "newton") KitBase.UnstructVSpace(-1, 1, 16, rand(16), ones(16)) diff --git a/test/test_theory.jl b/test/test_theory.jl index c21245d0..a73bb716 100644 --- a/test/test_theory.jl +++ b/test/test_theory.jl @@ -252,7 +252,7 @@ KitBase.boltzmann_fft!(rand(16, 16, 16), rand(16, 16, 16), fsm) KitBase.boltzmann_ode!(zeros(16, 16, 16), rand(16, 16, 16), (1.0, 5, phi, psi, phipsi), 0.0) KitBase.bgk_ode!(zeros(16, 16, 16), rand(16, 16, 16), (rand(16, 16, 16), 1e-2), 0.0) -vs = KitBase.VSpace3D(-5.0, 5.0, 16, -5.0, 5.0, 16, -5.0, 5.0, 16, "algebra") +vs = KitBase.VSpace3D(-5.0, 5.0, 16, -5.0, 5.0, 16, -5.0, 5.0, 16, type = "algebra") u, v, w = vs.u[:, 1, 1], vs.v[1, :, 1], vs.w[1, 1, :] vnu = hcat(vs.u[:], vs.v[:], vs.w[:]) uuni1d = linspace(vs.u[1, 1, 1], vs.u[end, 1, 1], 16) From 9d1115826aeb96ada902511aa5b8f879ec5bb547 Mon Sep 17 00:00:00 2001 From: vavrines Date: Fri, 5 Aug 2022 11:32:58 +0200 Subject: [PATCH 2/5] Add f_maxwellian methods --- src/Theory/theory_maxwellian.jl | 23 +++++++++++++++++++++++ test/test_theory.jl | 2 ++ 2 files changed, 25 insertions(+) diff --git a/src/Theory/theory_maxwellian.jl b/src/Theory/theory_maxwellian.jl index fb7e48bd..a219cd88 100644 --- a/src/Theory/theory_maxwellian.jl +++ b/src/Theory/theory_maxwellian.jl @@ -286,6 +286,7 @@ function mixture_maxwellian!(M::AM, u::T, v::T, w::T, prim::AM) where {T<:AM} end + """ $(SIGNATURES) @@ -308,3 +309,25 @@ function mixture_energy_maxwellian(h::AA{T,3}, prim::AM, K) where {T} return b end + + +""" +$(SIGNATURES) + +Reduced Maxwellian distribution related to energy for mixture +""" +function f_maxwellian(f::AV, u::AV, weights::AV, γ = 3) + w = moments_conserve(f, u, weights) + prim = conserve_prim(w, γ) + + return maxwellian(u, prim) +end + +f_maxwellian(f::AV, vs = VSpace1D(-6, 6, size(f, 1); precision = Float32)::VSpace1D, γ = 3) = + f_maxwellian(f, vs.u, vs.weights, γ) + +function f_maxwellian(f::AM, vs = VSpace1D(-6, 6, size(f, 1); precision = Float32)::VSpace1D, γ = 3) + M = [f_maxwellian(f[:, i], vs, γ) for i in axes(f, 2)] + + return hcat(M...) +end diff --git a/test/test_theory.jl b/test/test_theory.jl index a73bb716..f092753d 100644 --- a/test/test_theory.jl +++ b/test/test_theory.jl @@ -221,6 +221,8 @@ KitBase.rykov!( 0.3049, ) +KitBase.f_maxwellian(rand(16, 10)) + KitBase.reduce_distribution(randn(16, 51), ω, 1) KitBase.reduce_distribution(randn(16, 24, 24), ones(24, 24), 1) KitBase.reduce_distribution( From a85bcc66111d949dfab3a86a2b06b9c5d056a769 Mon Sep 17 00:00:00 2001 From: Tianbai Xiao Date: Fri, 5 Aug 2022 11:46:04 +0200 Subject: [PATCH 3/5] Export f_maxwellian --- src/Theory/theory.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Theory/theory.jl b/src/Theory/theory.jl index fccf8886..07a6ced8 100755 --- a/src/Theory/theory.jl +++ b/src/Theory/theory.jl @@ -9,7 +9,7 @@ export gauss_moments, mixture_gauss_moments, discrete_moments export moments_conserve, diatomic_moments_conserve, mixture_moments_conserve, flux_conserve! export pdf_slope, mixture_pdf_slope, moments_conserve_slope, mixture_moments_conserve_slope export pressure, stress, heat_flux -export maxwellian, energy_maxwellian, maxwellian! +export maxwellian, energy_maxwellian, maxwellian!, f_maxwellian export mixture_maxwellian, mixture_energy_maxwellian, mixture_maxwellian! export shakhov, shakhov!, rykov! export reduce_distribution, full_distribution, shift_pdf! From e0c49f5b6a53b1b797e8d0b2072052f3c6e07841 Mon Sep 17 00:00:00 2001 From: Tianbai Xiao Date: Fri, 5 Aug 2022 21:06:00 +0200 Subject: [PATCH 4/5] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 4a13fd84..0b6eeeab 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "KitBase" uuid = "86eed249-3a28-466f-8d3a-596821e1af9a" authors = ["Tianbai Xiao "] -version = "0.9.2" +version = "0.9.3" [deps] BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0" From 5b71751ec627ea936a316f819051b9014eb8561e Mon Sep 17 00:00:00 2001 From: vavrines Date: Fri, 5 Aug 2022 21:22:46 +0200 Subject: [PATCH 5/5] Clean shakhov codes --- src/Theory/theory_shakhov.jl | 78 ++++++++++++++++++++++-------------- 1 file changed, 48 insertions(+), 30 deletions(-) diff --git a/src/Theory/theory_shakhov.jl b/src/Theory/theory_shakhov.jl index 2b4cc92c..50256500 100644 --- a/src/Theory/theory_shakhov.jl +++ b/src/Theory/theory_shakhov.jl @@ -2,8 +2,10 @@ $(SIGNATURES) Shakhov non-equilibrium part + +1F1V """ -function shakhov(u::AV{X}, M::AV{Y}, q, prim::AV{Z}, Pr) where {X<:FN,Y<:FN,Z<:Real} # 1F1V +function shakhov(u::AV, M::AV, q, prim::AV, Pr) M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * (u - prim[2]) * q * @@ -13,16 +15,20 @@ function shakhov(u::AV{X}, M::AV{Y}, q, prim::AV{Z}, Pr) where {X<:FN,Y<:FN,Z<:R return M_plus end -#--- 2F1V ---# +""" +$(SIGNATURES) + +2F1V +""" function shakhov( - u::AV{X}, - H::Y, - B::Y, + u::AV, + H::T, + B::T, q, - prim::AV{Z}, + prim::AV, Pr, K, -) where {X<:FN,Y<:AV{<:FN},Z<:Real} +) where {T<:AV} H_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * (u - prim[2]) * @@ -39,15 +45,19 @@ function shakhov( end -#--- 1F2V ---# +""" +$(SIGNATURES) + +1F2V +""" function shakhov( u::T, v::T, - M::AM{X}, - q::AV{Y}, - prim::AV{Z}, + M::AM, + q::AV, + prim::AV, Pr, -) where {T<:AA{<:FN,2},X<:FN,Y<:FN,Z<:Real} +) where {T<:AM} M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) * @@ -58,17 +68,21 @@ function shakhov( end -#--- 2F2V ---# +""" +$(SIGNATURES) + +2F2V +""" function shakhov( u::T, v::T, H::X, B::X, - q::AV{Y}, - prim::AV{Z}, + q::AV, + prim::AV, Pr, K, -) where {T<:AA{<:FN,2},X<:AA{<:FN,2},Y<:Real,Z<:Real} +) where {T<:AM,X<:AM} H_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * ((u - prim[2]) * q[1] + (v - prim[3]) * q[2]) * @@ -83,16 +97,20 @@ function shakhov( end -#--- 1F3V ---# +""" +$(SIGNATURES) + +1F3V +""" function shakhov( u::T, v::T, w::T, - M::AA{X,3}, - q::AV{Y}, - prim::AV{Z}, + M::AA{T1,3}, + q::AV, + prim::AV, Pr, -) where {T<:AA{<:FN,3},X<:FN,Y<:Real,Z<:Real} +) where {T<:AA{T2,3},T1} where {T2} M_plus = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * ((u - prim[2]) * q[1] + (v - prim[3]) * q[2] + (w - prim[4]) * q[3]) * @@ -111,7 +129,7 @@ In-place Shakhov non-equilibrium part 1F1V """ -function shakhov!(S::T1, u::AV{T2}, M::T1, q, prim, Pr) where {T1<:AA{<:FN,1},T2<:FN} +function shakhov!(S::T1, u::AV, M::T1, q, prim, Pr) where {T1<:AV} @. S = @. 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * (u - prim[2]) * q * @@ -129,14 +147,14 @@ $(SIGNATURES) function shakhov!( SH::T1, SB::T1, - u::AV{T2}, + u::AV, H::T1, B::T1, q, prim, Pr, K, -) where {T1<:AA{<:FN,1},T2<:FN} +) where {T1<:AV} @. SH = 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * @@ -165,10 +183,10 @@ function shakhov!( u::T2, v::T2, M::T1, - q::AV{T3}, + q::AV, prim, Pr, -) where {T1<:AA{<:FN,2},T2<:AA{<:FN,2},T3<:Real} +) where {T1<:AM,T2<:AM} @. S = 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * @@ -192,11 +210,11 @@ function shakhov!( v::T2, H::T1, B::T1, - q::AV{T3}, + q::AV, prim, Pr, K, -) where {T1<:AA{<:FN,2},T2<:AA{<:FN,2},T3<:Real} +) where {T1<:AM,T2<:AM} @. SH = 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] * @@ -224,10 +242,10 @@ function shakhov!( v::T2, w::T2, M::T1, - q::AV{T3}, + q::AV, prim, Pr, -) where {T1<:AA{<:FN,3},T2<:AA{<:FN,3},T3<:Real} +) where {T1<:AA{T3,3},T2<:AA{T4,3}} where {T3,T4} @. S = 0.8 * (1.0 - Pr) * prim[end]^2 / prim[1] *