Anisotropy (#63)

* Add anisotropy computation

* Add LinearAlgebra, StaticArrays in deps

* Bump patch

Project.toml

@@ -1,21 +1,23 @@
 name = "Batsrus"
 uuid = "e74ebddf-6ac1-4047-a0e5-c32c99e57753"
 authors = ["Hongyang Zhou <[email protected]>"]
-version = "0.6.7"
+version = "0.6.8"
 
 [deps]
 FortranFiles = "c58ffaec-ab22-586d-bfc5-781a99fd0b10"
 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"
 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"

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
 

src/Batsrus.jl

@@ -3,10 +3,12 @@ module Batsrus
 #
 # Hongyang Zhou, [email protected]
 
+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")

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

src/plot/utility.jl → 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

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"
+    sha256 = "e344560deced105e684ed975276f34d7003427f5337c7e4b1def8274e5470a23"

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