Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

E calculation #71

Merged
merged 6 commits into from
Dec 6, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
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
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
Original file line number Diff line number Diff line change
Expand Up @@ -188,9 +188,14 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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

Check warning on line 314 in src/select.jl

View check run for this annotation

Codecov / codecov/patch

src/select.jl#L314

Added line #L314 was not covered by tests

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)),

Check warning on line 350 in src/select.jl

View check run for this annotation

Codecov / codecov/patch

src/select.jl#L349-L350

Added lines #L349 - L350 were not covered by tests
"B" => (bd -> fill_vector_from_scalars(bd, :B)),
"E" => (bd -> fill_vector_from_scalars(bd, :E)),
"U" => (bd -> fill_vector_from_scalars(bd, :U)),

Check warning on line 353 in src/select.jl

View check run for this annotation

Codecov / codecov/patch

src/select.jl#L353

Added line #L353 was not covered by tests
"Anisotropy0" => (bd -> get_anisotropy(bd, 0)),
"Anisotropy1" => (bd -> get_anisotropy(bd, 1)),
)
37 changes: 37 additions & 0 deletions src/utility.jl
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,43 @@
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)

Check warning on line 121 in src/utility.jl

View check run for this annotation

Codecov / codecov/patch

src/utility.jl#L120-L121

Added lines #L120 - L121 were not covered by tests
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
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
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
Loading