Improve tec conversion (#54)

* Update docs * Change tec convert API * Bump patch
henry2004y · Jun 10, 2024 · d5f49f3 · d5f49f3 · henry2004y · Jun 10, 2024
1 parent 77ebf59
commit d5f49f3
Show file tree

Hide file tree

Showing 5 changed files with 37 additions and 58 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Batsrus"
 uuid = "e74ebddf-6ac1-4047-a0e5-c32c99e57753"
 authors = ["Hongyang Zhou <[email protected]>"]
-version = "0.5.8"
+version = "0.5.9"
 
 [deps]
 FortranFiles = "c58ffaec-ab22-586d-bfc5-781a99fd0b10"

diff --git a/docs/src/man/manual.md b/docs/src/man/manual.md
@@ -59,22 +59,21 @@ Here is a full list of predefined derived quantities:
 
 We can convert 2D/3D BATSRUS outputs `*.dat` to VTK formats. It uses the VTK XML format writer [writeVTK](https://github.com/jipolanco/WriteVTK.jl) to generate files for Paraview and Tecplot. The default converted filename is `out.vtu`.
 
-ASCII Tecplot file (supports both `tec` and `tcp`) and binary Tecplot file (set `DOSAVETECBINARY=TRUE` in BATSRUS `PARAM.in`):
+- ASCII Tecplot file (supports both `tec` and `tcp`) and binary Tecplot file (set `DOSAVETECBINARY=TRUE` in BATSRUS `PARAM.in`):
 
 ```julia
 file = "x=0_mhd_1_n00000050.dat"
-head, data, connectivity = readtecdata(file)
-convertTECtoVTU(head, data, connectivity)
+convertTECtoVTU(file)
 ```
 
-3D structured IDL file (`gridType=1` returns rectilinear `vtr` file, `gridType=2` returns structured `vts` file):
+- 3D structured IDL file (`gridType=1` returns rectilinear `vtr` file, `gridType=2` returns structured `vts` file):
 
 ```julia
 file = "3d_structured.out"
 convertIDLtoVTK(file, gridType=1)
 ```
 
-3D unstructured IDL file together with header and tree file:
+- 3D unstructured IDL file together with header and tree file:
 
 ```julia
 filetag = "3d_var_1_n00002500"
@@ -84,52 +83,28 @@ convertIDLtoVTK(filetag)
 !!! note
     The file suffix should not be provided for this to work correctly!
 
-Multiple files:
+- Multiple files:
 
 ```julia
-using Batsrus, Glob
-filenamesIn = "3d*.dat"
 dir = "."
-filenames = Vector{String}(undef, 0)
-filesfound = glob(filenamesIn, dir)
-filenames = vcat(filenames, filesfound)
-tec = readtecdata.(filenames) # head, data, connectivity
-for (i, outname) in enumerate(filenames)
-   convertTECtoVTU(tec[i][1], tec[i][2], tec[i][3], outname[1:end-4])
-end
-```
+filenames = filter(file -> startswith(file, "3d") && endswith(file, ".dat"), readdir(dir))
+filenames = dir .* filenames
 
-If each individual file size is large, consider using:
-
-```julia
-using Batsrus, Glob
-filenamesIn = "3d*.dat"
-dir = "."
-filenames = Vector{String}(undef, 0)
-filesfound = glob(filenamesIn, dir)
-filenames = vcat(filenames, filesfound)
-for (i, outname) in enumerate(filenames)
-   head, data, connectivity = readtecdata(outname)
-   convertTECtoVTU(head, data, connectivity, outname[1:end-4])
+for filename in filenames
+   convertTECtoVTU(filename, filename[1:end-4])
 end
 ```
 
-Multiple files in parallel:
+* Processing multiple files with threads in parallel:
 
 ```julia
-using Distributed
-@everywhere using Batsrus, Glob
-
-filenamesIn = "cut*.dat"
 dir = "."
-filenames = Vector{String}(undef, 0)
-filesfound = glob(filenamesIn, dir)
-filenames = vcat(filenames, filesfound)
-
-@sync @distributed for outname in filenames
-   println("filename=$(outname)")
-   head, data, connectivity = readtecdata(outname)
-   convertTECtoVTU(head, data, connectivity, outname[1:end-4])
+filenames = filter(file -> startswith(file, "3d") && endswith(file, ".dat"), readdir(dir))
+filenames = dir .* filenames
+
+Threads.@threads for filename in filenames
+   println("filename=$filename")
+   convertTECtoVTU(filename, filename[1:end-4])
 end
 ```
 

diff --git a/src/io.jl b/src/io.jl
@@ -62,8 +62,8 @@ function readlogdata(file::AbstractString)
    nx        = 1
    nw        = length(variables)
 
-   data = zeros(nw,nLine)
-   for i in 1:nLine
+   data = zeros(nw, nLine)
+   @inbounds for i in 1:nLine
       line = split(readline(f))
       data[:,i] = parse.(Float64, line)
    end
@@ -411,7 +411,7 @@ end
 function getascii!(x, w, fileID::IOStream, filehead::NamedTuple)
    ndim = filehead.ndim
    Ids = CartesianIndices(size(x)[1:ndim])
-   for ids in Ids
+   @inbounds @views for ids in Ids
       temp = parse.(Float64, split(readline(fileID)))
       x[ids,:] = temp[1:ndim]
       w[ids,:] = temp[ndim+1:end]
@@ -425,7 +425,7 @@ function getbinary!(x, w, fileID::IOStream, filehead::NamedTuple)
    dimlast = filehead.ndim + 1
    read!(fileID, x)
    skip(fileID, 2*TAG)
-   for iw in axes(w, dimlast)
+   @inbounds for iw in axes(w, dimlast)
       read!(fileID, selectdim(w, dimlast, iw))
       skip(fileID, 2*TAG)
    end

diff --git a/src/vtk.jl b/src/vtk.jl
@@ -69,16 +69,17 @@ struct Head
    dxPlot_D::Vector{Float64}
 end
 
+"BATSRUS output high-level struct."
 struct Batl
    head::Head
    iTree_IA::Array{Int32,2}
    iRatio_D::Vector{Int32} # Array of refinement ratios
    nDim::Int8
 end
 
-
 """
-    convertTECtoVTU(head, data, connectivity, filename="out")
+    convertTECtoVTU(file::AbstractString, outname="out")
+    convertTECtoVTU(head, data, connectivity, outname="out")
 
 Convert unstructured Tecplot data to VTK. Note that if using voxel type data in VTK, the
 connectivity sequence is different from Tecplot: the 3D connectivity sequence in Tecplot is
@@ -92,7 +93,12 @@ for i in 1:2
 end
 ```
 """
-function convertTECtoVTU(head, data, connectivity, filename="out")
+function convertTECtoVTU(file::AbstractString, outname="out")
+   head, data, connectivity = readtecdata(file)
+   convertTECtoVTU(head, data, connectivity, outname)
+end
+
+function convertTECtoVTU(head, data, connectivity, outname="out")
    nVar = length(head.variables)
    points = @view data[1:head.nDim,:]
    cells = Vector{MeshCell{VTKCellType,Array{Int32,1}}}(undef, head.nCell)
@@ -107,7 +113,7 @@ function convertTECtoVTU(head, data, connectivity, filename="out")
       end
    end
 
-   vtkfile = vtk_grid(filename, points, cells)
+   vtkfile = vtk_grid(outname, points, cells)
 
    for ivar in head.nDim+1:nVar
       if endswith(head.variables[ivar],"_x") # vector
@@ -123,7 +129,7 @@ function convertTECtoVTU(head, data, connectivity, filename="out")
             namevar = replace(head.variables[ivar], "_x"=>"")
             vtk_point_data(vtkfile, (var1, var2), namevar)
          end
-      elseif endswith(head.variables[ivar],r"_y|_z")
+      elseif endswith(head.variables[ivar], r"_y|_z")
          continue
       else
          var = @view data[ivar,:]
@@ -161,7 +167,7 @@ function convertIDLtoVTK(filename::AbstractString; gridType::Int=1, verbose::Boo
          z = @view data.x[1,1,:,3]
 
          outfiles = vtk_grid(outname, x,y,z) do vtk
-            for ivar = 1:nVar
+            for ivar in 1:nVar
                if data.head.wnames[ivar][end] == 'x' # vector
                   var1 = @view data.w[:,:,:,ivar]
                   var2 = @view data.w[:,:,:,ivar+1]
@@ -180,7 +186,7 @@ function convertIDLtoVTK(filename::AbstractString; gridType::Int=1, verbose::Boo
          xyz = permutedims(data.x, [4,1,2,3])
 
          outfiles = vtk_grid(outname, xyz) do vtk
-            for ivar = 1:nVar
+            for ivar in 1:nVar
                if data.head.wnames[ivar][end] == 'x' # vector
                   var1 = @view data.w[:,:,:,ivar]
                   var2 = @view data.w[:,:,:,ivar+1]
@@ -732,7 +738,7 @@ function getConnectivity(batl::Batl)
    # Pre-allocate just to let Julia know this variable.
    connectivity = Array{Int32,2}(undef, nConn, nElem)
 
-   for iRound in 1:2
+   @inbounds for iRound in 1:2
       if iRound == 2
          connectivity = Array{Int32,2}(undef, nConn, nElem)
          iElem = 0
@@ -751,7 +757,7 @@ function getConnectivity(batl::Batl)
             # initial cell index
             iCell = nI*nJ*nK*(nBlockBefore + iBlock - 1)
 
-            @inbounds for k in 1:nK, j in 1:nJ, i in 1:nI
+            for k in 1:nK, j in 1:nJ, i in 1:nI
                iCell += 1
                iCell_G[i+1,j+1,k+1] = iCell
             end

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -82,9 +82,7 @@ end
 
    @testset "VTK" begin
       file = joinpath(datapath, "3d_bin.dat")
-      head, bd, connectivity = readtecdata(file)
-      @test maximum(connectivity) ≤ head[:nNode] # check if it's read correctly
-      convertTECtoVTU(head, bd, connectivity)
+      convertTECtoVTU(file)
       sha_str = bytes2hex(open(sha1, "out.vtu"))
       @test sha_str == "5b04747666542d802357dec183177f757754a254"
       rm("out.vtu")