Skip to content

Commit

Permalink
E calculation (#71)
Browse files Browse the repository at this point in the history 
* Add E field computations; Add axis extraction

* Export more utility methods

* Bump patch
  • Loading branch information
@henry2004y
henry2004y authored Dec 6, 2024
1 parent 62897ff commit b11ecb3
Show file tree
Hide file tree
Showing 6 changed files with 141 additions and 28 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Batsrus"
uuid = "e74ebddf-6ac1-4047-a0e5-c32c99e57753"
authors = ["Hongyang Zhou <[email protected]>"]
version = "0.7.2"
version = "0.7.3"

[deps]
FortranFiles = "c58ffaec-ab22-586d-bfc5-781a99fd0b10"
Expand Down
4 changes: 3 additions & 1 deletion src/Batsrus.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
114 changes: 89 additions & 25 deletions src/select.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)),
)
"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)),
)
37 changes: 37 additions & 0 deletions src/utility.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion src/vtk.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
10 changes: 10 additions & 0 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,19 +69,29 @@ 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]
@test bd["Ue2"][2,1] == 33784.973f0
@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
Expand Down

2 comments on commit b11ecb3

@henry2004y
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/120870

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.7.3 -m "<description of version>" b11ecb3efb451f568237ca8bb61986c6f07eebe6
git push origin v0.7.3

Please sign in to comment.