diff --git a/Project.toml b/Project.toml index c48bd5d3..0a49aa99 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Batsrus" uuid = "e74ebddf-6ac1-4047-a0e5-c32c99e57753" authors = ["Hongyang Zhou "] -version = "0.6.7" +version = "0.6.8" [deps] FortranFiles = "c58ffaec-ab22-586d-bfc5-781a99fd0b10" @@ -9,6 +9,7 @@ Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NaturalNeighbours = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6" Parsers = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -16,6 +17,7 @@ PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" @@ -38,6 +40,7 @@ PyPlot = "2.9" RecipesBase = "1.1" Reexport = "1.2" Requires = "1.1" +StaticArrays = "1.2" Unitful = "1.7" WriteVTK = "1.9" julia = "1.6" diff --git a/docs/src/man/manual.md b/docs/src/man/manual.md index c71cc945..4d02b3e7 100644 --- a/docs/src/man/manual.md +++ b/docs/src/man/manual.md @@ -66,15 +66,16 @@ Here is a full list of predefined derived quantities: | Derived variable name | Meaning | Required variable | |-----------------------|----------------------------------|-------------------| -| B2 | magnetic field magnitude squared | Bx, By, Bz | -| E2 | electric field magnitude squared | Ex, Ey, Ez | -| U2 | velocity magnitude squared | Ux, Uy, Uz | -| Bmag | magnetic field magnitude | Bx, By, Bz | -| Emag | electric field magnitude | Ex, Ey, Ez | -| Umag | velocity magnitude | Ux, Uy, Uz | -| B | magnetic field vector | Bx, By, Bz | -| E | electric field vector | Ex, Ey, Ez | -| U | velocity vector | Ux, Uy, Uz | +| B2 | Magnetic field magnitude squared | Bx, By, Bz | +| E2 | Electric field magnitude squared | Ex, Ey, Ez | +| U2 | Velocity magnitude squared | Ux, Uy, Uz | +| Bmag | Magnetic field magnitude | Bx, By, Bz | +| Emag | Electric field magnitude | Ex, Ey, Ez | +| Umag | Velocity magnitude | Ux, Uy, Uz | +| B | Magnetic field vector | Bx, By, Bz | +| E | Electric field vector | Ex, Ey, Ez | +| U | Velocity vector | Ux, Uy, Uz | +| Anisotropy | Pressure/Temperature anisotropy | B, P tensor | ## Output format conversion diff --git a/src/Batsrus.jl b/src/Batsrus.jl index 8b5fa5f0..51517f68 100644 --- a/src/Batsrus.jl +++ b/src/Batsrus.jl @@ -3,10 +3,12 @@ module Batsrus # # Hongyang Zhou, hyzhou@umich.edu +using LinearAlgebra: normalize, ×, ⋅ using Printf, Reexport, Requires using Parsers using Interpolations: cubic_spline_interpolation, BSpline, Linear, scale, interpolate import NaturalNeighbours as NN +using StaticArrays: SVector, @SMatrix, SA export BATLData, load, readlogdata, readtecdata, showhead, # io @@ -48,7 +50,7 @@ using .UnitfulBatsrus include("io.jl") include("select.jl") include("vtk.jl") -include("plot/utility.jl") +include("utility.jl") include("plot/plots.jl") include("hdf.jl") diff --git a/src/select.jl b/src/select.jl index 28f7a7d0..41cae86e 100644 --- a/src/select.jl +++ b/src/select.jl @@ -221,6 +221,33 @@ function _fill_vector_from_scalars(bd::BATLData{ndim, T}, vstr1, vstr2, vstr3) w v end +function _get_anisotropy(bd::BATLData{2, T}, species=0) where {T} + Bx, By, Bz = bd["Bx"], bd["By"], bd["Bz"] + # Rotate the pressure tensor to align the 3rd direction with B + pop = string(species) + Pxx = bd["pXXS"*pop] + Pyy = bd["pYYS"*pop] + Pzz = bd["pZZS"*pop] + Pxy = bd["pXYS"*pop] + Pxz = bd["pXZS"*pop] + Pyz = bd["pYZS"*pop] + + Paniso = similar(Pxx) + + @inbounds for j in axes(Pxx, 2), i in axes(Pxx, 1) + b̂ = normalize(SA[Bx[i,j], By[i,j], Bz[i,j]]) + P = @SMatrix [ + Pxx[i,j] Pxy[i,j] Pxz[i,j]; + Pxy[i,j] Pyy[i,j] Pyz[i,j]; + Pxz[i,j] Pyz[i,j] Pzz[i,j]] + + Prot = rotateTensorToVectorZ(P, b̂) + Paniso[i,j] = (Prot[1,1] + Prot[2,2]) / (2*Prot[3,3]) + end + + Paniso +end + # Define derived parameters const variables_predefined = Dict( "B2" => (bd -> @. $getvar(bd, "Bx")^2 + $getvar(bd, "By")^2 + $getvar(bd, "Bz")^2), @@ -234,4 +261,6 @@ const variables_predefined = Dict( "B" => (bd -> _fill_vector_from_scalars(bd, "Bx", "By", "Bz")), "E" => (bd -> _fill_vector_from_scalars(bd, "Ex", "Ey", "Ez")), "U" => (bd -> _fill_vector_from_scalars(bd, "Ux", "Uy", "Uz")), + "Anisotropy0" => (bd -> _get_anisotropy(bd, 0)), + "Anisotropy1" => (bd -> _get_anisotropy(bd, 1)), ) \ No newline at end of file diff --git a/src/plot/utility.jl b/src/utility.jl similarity index 82% rename from src/plot/utility.jl rename to src/utility.jl index 3af8311d..9f37e42b 100644 --- a/src/plot/utility.jl +++ b/src/utility.jl @@ -1,4 +1,4 @@ -# Utility functions for plotting. +# Utility functions for plotting and analyzing. """ interp2d(bd::BATLData, var::AbstractString, plotrange=[-Inf, Inf, -Inf, Inf], @@ -205,4 +205,47 @@ function get_range(bd::BATLData, var::AbstractString, verbose=true) varIndex_ = findindex(bd, var) w = selectdim(bd.w, bd.head.ndim+1, varIndex_) wmin, wmax = extrema(w) +end + +""" + rotateTensorToVectorZ(tensor::AbstractMatrix, v::AbstractVector) -> SMatrix{3,3} + +Rotate `tensor` with a rotation matrix that aligns the 3rd direction with `vector`, which is equivalent to change the basis from (i,j,k) to (i′,j′,k′) where k′ ∥ vector. +Reference: [Tensor rotation](https://math.stackexchange.com/questions/2303869/tensor-rotation) +""" +function rotateTensorToVectorZ(tensor::AbstractMatrix{T}, v) where T + k = SVector{3, T}(0.0, 0.0, 1.0) + axis = v × k |> normalize + if axis[1] == axis[2] == 0 + return tensor + else + angle = acos(v ⋅ k / sqrt(v[1]^2 + v[2]^2 + v[3]^2)) + R = getRotationMatrix(axis, angle) + return R * tensor * R' + end +end + +""" + getRotationMatrix(axis::AbstractVector, angle::Real) --> SMatrix{3,3} + +Create a rotation matrix for rotating a 3D vector around a unit `axis` by an `angle` in radians. +Reference: [Rotation matrix from axis and angle](https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle) + +# Example + +```julia +using LinearAlgebra +v = [-0.5, 1.0, 1.0] +v̂ = normalize(v) +θ = deg2rad(-74) +R = getRotationMatrix(v̂, θ) +``` +""" +function getRotationMatrix(v::AbstractVector{<:AbstractFloat}, θ::Real) + sinθ, cosθ = sincos(eltype(v)(θ)) + tmp = 1 - cosθ + m = @SMatrix [ + cosθ+v[1]^2*tmp v[1]*v[2]*tmp-v[3]*sinθ v[1]*v[3]*tmp+v[2]*sinθ; + v[1]*v[2]*tmp+v[3]*sinθ cosθ+v[2]^2*tmp v[2]*v[3]*tmp-v[1]*sinθ; + v[1]*v[3]*tmp-v[2]*sinθ v[3]*v[2]*tmp+v[1]*sinθ cosθ+v[3]^2*tmp] end \ No newline at end of file diff --git a/test/Artifacts.toml b/test/Artifacts.toml index b3206126..51f0e5e3 100644 --- a/test/Artifacts.toml +++ b/test/Artifacts.toml @@ -1,7 +1,7 @@ [testdata] -git-tree-sha1 = "ecff1c52f54b1847036f98d8000b7268a5c3c935" +git-tree-sha1 = "d884abbeac09834475c56b256659becc1da64927" lazy = true [[testdata.download]] url = "https://github.com/henry2004y/batsrus_data/raw/master/batsrus_data.tar.gz" - sha256 = "fb2992dd3bca73162086be79c763603139740408ce87bbfe478593356d270600" \ No newline at end of file + sha256 = "e344560deced105e684ed975276f34d7003427f5337c7e4b1def8274e5470a23" \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 49ec762d..7b4d841f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -64,6 +64,11 @@ end w = slice1d(bd, "rho", 1, 1) @test sum(w) == 4.0f0 @test get_var_range(bd, "rho") == (0.11626893f0, 1.0f0) + + file = "z=0_fluid_region0_0_t00001640_n00010142.out" + bd = load(joinpath(datapath, file)) + @test bd["Anisotropy0"][1:2,1] ≈ Float32[1.2630985, 2.4700143] + @test bd["Anisotropy1"][1:2,1] ≈ Float32[1.2906302, 2.6070855] end @testset "Reading 2D unstructured" begin