diff --git a/Project.toml b/Project.toml index a75d8ac9..4ed6c994 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Batsrus" uuid = "e74ebddf-6ac1-4047-a0e5-c32c99e57753" authors = ["Hongyang Zhou "] -version = "0.7.2" +version = "0.7.3" [deps] FortranFiles = "c58ffaec-ab22-586d-bfc5-781a99fd0b10" diff --git a/src/Batsrus.jl b/src/Batsrus.jl index 731b8e3b..b56284a7 100644 --- a/src/Batsrus.jl +++ b/src/Batsrus.jl @@ -12,7 +12,9 @@ using StaticArrays: SVector, @SMatrix, SA export BATS, load, readlogdata, readtecdata, showhead, # io - getvar, cutdata, subvolume, subsurface, # select + getvar, cutdata, subvolume, subsurface, get_convection_E, get_hall_E, + get_anisotropy, get_vectors, get_magnitude, get_magnitude2, + fill_vector_from_scalars, # select Batl, convertTECtoVTU, convertIDLtoVTK, readhead, readtree, getConnectivity, # vtk interp1d, interp2d, slice1d, get_var_range, squeeze, get_range # plot/utility diff --git a/src/select.jl b/src/select.jl index a4ab4951..106a7211 100644 --- a/src/select.jl +++ b/src/select.jl @@ -188,9 +188,14 @@ end getvar(bd, var) -"Construct vectors from scalar components." -function _fill_vector_from_scalars(bd::BATS{ndim, ndimp1, T}, var) where {ndim, ndimp1, T} - v1, v2, v3 = _get_vectors(bd, var) +""" + fill_vector_from_scalars(bd::BATS, var) + +Construct vector of `var` from its scalar components. +Alternatively, use [`get_vectors`](@ref) for returning the individual components. +""" +function fill_vector_from_scalars(bd::BATS{ndim, ndimp1, T}, var) where {ndim, ndimp1, T} + v1, v2, v3 = get_vectors(bd, var) v = Array{T, ndimp1}(undef, 3, size(v1)...) Rpost = CartesianIndices(size(v1)) @@ -203,8 +208,13 @@ function _fill_vector_from_scalars(bd::BATS{ndim, ndimp1, T}, var) where {ndim, v end -function _get_magnitude2(bd::BATS{2, 3, T}, var=:B) where T - vx, vy, vz = _get_vectors(bd, var) +""" + get_magnitude2(bd::BATS, var) + +Calculate the magnitude square of vector `var`. See [`get_vectors`](@ref) for the options. +""" +function get_magnitude2(bd::BATS{2, 3, T}, var=:B) where T + vx, vy, vz = get_vectors(bd, var) v = similar(vx)::Array{T, 2} @inbounds @simd for i in eachindex(v) @@ -214,8 +224,13 @@ function _get_magnitude2(bd::BATS{2, 3, T}, var=:B) where T v end -function _get_magnitude(bd::BATS{2, 3, T}, var=:B) where T - vx, vy, vz = _get_vectors(bd, var) +""" + get_magnitude(bd::BATS, var) + +Calculate the magnitude of vector `var`. See [`get_vectors`](@ref) for the options. +""" +function get_magnitude(bd::BATS{2, 3, T}, var=:B) where T + vx, vy, vz = get_vectors(bd, var) v = similar(vx)::Array{T, 2} @inbounds @simd for i in eachindex(v) @@ -225,7 +240,12 @@ function _get_magnitude(bd::BATS{2, 3, T}, var=:B) where T v end -function _get_vectors(bd::BATS{Dim, Dimp1, T}, var) where {Dim, Dimp1,T} +""" + get_vectors(bd::BATS, var) + +Return a tuple of vectors of `var`. `var` can be `:B`, `:E`, `:U`, `:U0`, or `:U1`. +""" +function get_vectors(bd::BATS{Dim, Dimp1, T}, var) where {Dim, Dimp1, T} VT = SubArray{T, Dim, Array{T, Dimp1}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, true} if var == :B @@ -243,7 +263,12 @@ function _get_vectors(bd::BATS{Dim, Dimp1, T}, var) where {Dim, Dimp1,T} vx, vy, vz end -function _get_anisotropy(bd::BATS{2, 3, T}, species=0) where T +""" + get_anisotropy(bd::BATS, species=0) + +Calculate the pressure anisotropy for `species`, indexing from 0. +""" +function get_anisotropy(bd::BATS{2, 3, T}, species=0) where T VT = SubArray{T, 2, Array{T, 3}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, true} Bx, By, Bz = bd["Bx"]::VT, bd["By"]::VT, bd["Bz"]::VT @@ -272,21 +297,60 @@ function _get_anisotropy(bd::BATS{2, 3, T}, species=0) where T Paniso end +"Return the convection electric field from PIC outputs." +function get_convection_E(bd::BATS) + Bx, By, Bz = get_vectors(bd, :B) + # Let us use H+ velocities as the ion bulk velocity and ignore heavy ions + uix, uiy, uiz = get_vectors(bd, :U1) + + Econvx = similar(Bx) + Econvy = similar(By) + Econvz = similar(Bz) + # -Ui × B + @simd for i in eachindex(Econvx) + Econvx[i] = -uiy[i]*Bz[i] + uiz[i]*By[i] + Econvy[i] = -uiz[i]*Bx[i] + uix[i]*Bz[i] + Econvz[i] = -uix[i]*By[i] + uiy[i]*Bx[i] + end + + Econvx, Econvy, Econvz +end + +"Return the Hall electric field from PIC outputs." +function get_hall_E(bd::BATS) + Bx, By, Bz = get_vectors(bd, :B) + uex, uey, uez = get_vectors(bd, :U0) + # Let us use H+ velocities as the ion bulk velocity and ignore heavy ions + uix, uiy, uiz = get_vectors(bd, :U1) + + Ehallx = similar(Bx) + Ehally = similar(By) + Ehallz = similar(Bz) + # (Ui - Ue) × B + for i in eachindex(Ehallx) + Ehallx[i] = (uiy[i] - uey[i])*Bz[i] - (uiz[i] - uez[i])*By[i] + Ehally[i] = (uiz[i] - uez[i])*Bx[i] - (uix[i] - uex[i])*Bz[i] + Ehallz[i] = (uix[i] - uex[i])*By[i] - (uiy[i] - uey[i])*Bx[i] + end + + Ehallx, Ehally, Ehallz +end + # Define derived parameters const variables_predefined = Dict{String, Function}( - "B2" => (bd -> _get_magnitude2(bd, :B)), - "E2" => (bd -> _get_magnitude2(bd, :E)), - "U2" => (bd -> _get_magnitude2(bd, :U)), - "Ue2" => (bd -> _get_magnitude2(bd, :U0)), - "Ui2" => (bd -> _get_magnitude2(bd, :U1)), - "Bmag" => (bd -> _get_magnitude(bd, :B)), - "Emag" => (bd -> _get_magnitude(bd, :E)), - "Umag" => (bd -> _get_magnitude(bd, :U)), - "Uemag" => (bd -> _get_magnitude(bd, :U0)), - "Uimag" => (bd -> _get_magnitude(bd, :U1)), - "B" => (bd -> _fill_vector_from_scalars(bd, :B)), - "E" => (bd -> _fill_vector_from_scalars(bd, :E)), - "U" => (bd -> _fill_vector_from_scalars(bd, :U)), - "Anisotropy0" => (bd -> _get_anisotropy(bd, 0)), - "Anisotropy1" => (bd -> _get_anisotropy(bd, 1)), -) \ No newline at end of file + "B2" => (bd -> get_magnitude2(bd, :B)), + "E2" => (bd -> get_magnitude2(bd, :E)), + "U2" => (bd -> get_magnitude2(bd, :U)), + "Ue2" => (bd -> get_magnitude2(bd, :U0)), + "Ui2" => (bd -> get_magnitude2(bd, :U1)), + "Bmag" => (bd -> get_magnitude(bd, :B)), + "Emag" => (bd -> get_magnitude(bd, :E)), + "Umag" => (bd -> get_magnitude(bd, :U)), + "Uemag" => (bd -> get_magnitude(bd, :U0)), + "Uimag" => (bd -> get_magnitude(bd, :U1)), + "B" => (bd -> fill_vector_from_scalars(bd, :B)), + "E" => (bd -> fill_vector_from_scalars(bd, :E)), + "U" => (bd -> fill_vector_from_scalars(bd, :U)), + "Anisotropy0" => (bd -> get_anisotropy(bd, 0)), + "Anisotropy1" => (bd -> get_anisotropy(bd, 1)), +) diff --git a/src/utility.jl b/src/utility.jl index d36adaa0..3ff4bc55 100644 --- a/src/utility.jl +++ b/src/utility.jl @@ -89,6 +89,43 @@ function interp2d(bd::BATS{2, 3, T}, var::AbstractString, xi, yi, Wi end +"Return the axis range for 2D outputs. See [`interp2d`](@ref)." +function meshgrid(bd::BATS, + plotrange::Vector=[-Inf32, Inf32, -Inf32, Inf32], plotinterval::Real=Inf32) + x = bd.x + + if bd.head.gencoord # Generalized coordinates + X, Y = eachslice(x, dims=3) + X, Y = vec(X), vec(Y) + + adjust_plotrange!(plotrange, extrema(X), extrema(Y)) + # Set a heuristic value if not set + if isinf(plotinterval) + plotinterval = (plotrange[2] - plotrange[1]) / size(X, 1) + end + xi = range(plotrange[1], stop=plotrange[2], step=plotinterval) + yi = range(plotrange[3], stop=plotrange[4], step=plotinterval) + else # Cartesian coordinates + xrange = range(x[1,1,1], x[end,1,1], length=size(x,1)) + yrange = range(x[1,1,2], x[1,end,2], length=size(x,2)) + if all(isinf.(plotrange)) + xi, yi = xrange, yrange + else + adjust_plotrange!(plotrange, (xrange[1], xrange[end]), (yrange[1], yrange[end])) + + if isinf(plotinterval) + xi = range(plotrange[1], stop=plotrange[2], step=xrange[2] - xrange[1]) + yi = range(plotrange[3], stop=plotrange[4], step=yrange[2] - yrange[1]) + else + xi = range(plotrange[1], stop=plotrange[2], step=plotinterval) + yi = range(plotrange[3], stop=plotrange[4], step=plotinterval) + end + end + end + + xi, yi +end + "Find variable index in the BATSRUS data." function findindex(bd::BATS, var::AbstractString) varIndex_ = findfirst(x->lowercase(x)==lowercase(var), bd.head.wnames) diff --git a/src/vtk.jl b/src/vtk.jl index 05790c31..36931ff8 100644 --- a/src/vtk.jl +++ b/src/vtk.jl @@ -290,7 +290,7 @@ function swaprows!(X::Matrix, i::Int, j::Int) end end -"Return header from info file. Currently only designed for 2D and 3D." +"Return BATL header from info file. Currently only designed for 2D and 3D." function readhead(filehead) nDim = 3 nI, nJ, nK = 1, 1, 1 diff --git a/test/runtests.jl b/test/runtests.jl index 3c9d720e..4f8b95aa 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -69,6 +69,10 @@ end file = "z=0_fluid_region0_0_t00001640_n00010142.out" bd = load(joinpath(datapath, file)) + x, y = Batsrus.meshgrid(bd) + @test length(x) == 601 && y[2] == 0.0f0 + x, y = Batsrus.meshgrid(bd, Float32[-100, 100, -Inf, Inf]) + @test length(x) == 4 @test bd["Emag"][2,1] == 2655.4805f0 @test bd["E2"][2,1] == 7.051577f6 @test bd["E"][:,2,1] == Float32[-241.05942, -2644.2058, -40.53219] @@ -76,12 +80,18 @@ end @test bd["Ui2"][2,1] == bd["Ui2"][2,1] @test bd["Anisotropy0"][1:2,1] ≈ Float32[1.2630985, 2.4700143] @test bd["Anisotropy1"][1:2,1] ≈ Float32[1.2906302, 2.6070855] + w = get_convection_E(bd) + @test w[2][2,1] ≈ -2454.3933f0 + w = get_hall_E(bd) + @test w[2][2,1] ≈ -782.2945f0 end @testset "Reading 2D unstructured ascii" begin file = "bx0_mhd_6_t00000100_n00000352.out" bd = load(joinpath(datapath, file)) plotrange = [-Inf, Inf, -Inf, Inf] + x, y = Batsrus.meshgrid(bd) + @test length(x) == 117 && length(y) == 246 x, y, w = interp2d(bd, "rho", plotrange, useMatplotlib=false) @test w[1,2] == 5.000018304080387 @test bd["Umag"][2] == 71.85452748407637