# Help me optimize my special solution to save VTP (VTK) file?

**URL:** <https://discourse.julialang.org/t/help-me-optimize-my-special-solution-to-save-vtp-vtk-file/110064>\
**Category:** Performance\
**Created:** [February 11, 2024, 1:47pm UTC](https://discourse.julialang.org/t/help-me-optimize-my-special-solution-to-save-vtp-vtk-file/110064 "2024-02-11T13:47:15Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [February 11, 2024, 1:47pm UTC](https://discourse.julialang.org/t/help-me-optimize-my-special-solution-to-save-vtp-vtk-file/110064/1 "2024-02-11T13:47:15Z")

</div>

Hello!

Lately I’ve learned a lot from all of you guys on the forum on how to write performant code and reduce allocs/type instabilities, which has been really nice. I have a case where I want to save a lot of simulation date each N time steps, which is why I need my special solution instead of `WriteVTK` which is an absolutely beautiful library.

The code I have right now as testing is included here:

> **Summary**
>
> Current code, click to open
> 
> ```julia
> #https://www.analytech-solutions.com/analytech-solutions/blog/binary-io.html
> using XML
> using XML: Document, Declaration, Element, Text
> using StaticArrays
> 
> ### Functions=================================================
> # Function to create a DataArray element for VTK files
> function create_data_array_element(name::String, data::AbstractVector{T}, offset::Int) where T
> # Create the DataArray elements
> dataarray = Element("DataArray")
>     
> # Set attributes based on the input vector's type
> dataarray.attributes["type"] = string(eltype(first(data)))
> dataarray.attributes["Name"] = name  
> dataarray.attributes["NumberOfComponents"] = string(Int(sizeof(first(data))/sizeof(eltype(first(data)))))
> dataarray.attributes["format"] = "appended"
> dataarray.attributes["offset"] = string(offset)
> return dataarray
> end
> 
> ###===========================================================
> # Points = [SVector{3,Float64}(1,2,3), SVector{3,Float64}(4,5,6)]
> # Kernel = Float64.([100, 200]) 
> # KernelGradient = [SVector{3,Float64}(-1,1,0), SVector{3,Float64}(1,-1,0)]
> N = 6195
> Points = rand(SVector{3,Float64},N) * 10
> Kernel = rand(Float64,N) * 1000
> KernelGradient = rand(SVector{3,Float64},N) * 100
> 
> function PolyDataTemplate(filename::String, points, variable_names, args...)
> # Generate the XML document and then put in some fixed values
> xml_doc = Document(Declaration(version=1.0,encoding="utf-8"))
> vtk_file = Element("VTKFile")
> vtk_file.attributes["type"] = "PolyData"
> vtk_file.attributes["version"] = "1.0"
> vtk_file.attributes["byte_order"] = "LittleEndian"
> vtk_file.attributes["header_type"] = "UInt64"
> 
> # PolyData is the main section, filling it out
> polydata = Element("PolyData")
> piece = Element("Piece")
> N = length(points)
> piece.attributes["NumberOfPoints"] = string(N)
> 
> # This Points element and its associated DataArray has to be constructed individually
> points_element = Element("Points")
> point_dataarray = create_data_array_element("Points",points,0)
> point_dataarray["offset"] = 0
> 
>         
> # Generate appended data element
> appendeddata = Element("AppendedData")
> appendeddata.attributes["encoding"] = "raw"
> 
> # Start writing the file and generating the correct dataarrays with the right offsets in the loop
> NB = 0
> io = IOBuffer()
> write(io,"\n_")
> UncompressedHeaderN = N * length(first(points)) * sizeof(typeof(first(points)))
> NB += write(io, UncompressedHeaderN)
> NB += write(io, points)
> 
> # Generate XML tags for kwargs data
> pointdata = Element("PointData")
> dataarrays = Vector{XML.Node}(undef,length(args))
> 
> for i in eachindex(args)
> arg = args[i]
> dataarrays[i] = create_data_array_element(variable_names[i],arg,NB)
> 
> A = typeof(first(arg))
> T = eltype(A)
> Ni = length(arg)
> Tsz = sizeof(T)
> Nc = Int( sizeof(A) / Tsz )
> HowManyBytes = Tsz*Nc*Ni + Tsz
> 
> NB += HowManyBytes
> 
> write(io, NB)
> write(io, arg)
> end
> 
> # Take the result from the buffer, turn to string and write it
> v = take!(io)
> t = Text(String(v))
> write(io,"\n")
> push!(appendeddata,t)
> close(io)
> 
> # Glue all xml pieces together
> push!(xml_doc,vtk_file)
> push!(points_element,point_dataarray)
> push!(piece,points_element)
> push!(polydata,piece)
> push!(vtk_file,polydata)
> map(x -> push!(pointdata,x), dataarrays)
> push!(piece,pointdata)
> push!(vtk_file,appendeddata)
> 
> XML.write(filename,xml_doc)
> end
> 
> save_location = raw"E:\SPH\TestOfFile.vtp"
> 
> d = @report_opt target_modules=(@ __MODULE__ ,) PolyDataTemplate(save_location, Points, ["Kernel", "KernelGradient"], Kernel, KernelGradient)
> println(d)
> 
> @profview PolyDataTemplate(save_location, Points, ["Kernel", "KernelGradient"], Kernel, KernelGradient)
> 
> b = @benchmark PolyDataTemplate($save_location, $Points, $(["Kernel", "KernelGradient"]), $Kernel, $KernelGradient)
> display(b)
> 
> @code_warntype PolyDataTemplate(save_location, Points, ["Kernel", "KernelGradient"], Kernel, KernelGradient)
> 
> PolyDataTemplate(save_location, Points, ["Kernel", "KernelGradient"], Kernel, KernelGradient)
> 
> ```

It outputs the following which can be visualized in ParaView:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/3/131d7a0a4e2c3afde9f678ebcfed66fad665016b.png)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/0/f/0fc09559f2084f82a622633ebcfec6ab6c8d1b49.png)

I have tested my solution up against `WriteVTK` in some production code and found that my CustomVTP is indeed faster:

![image](https://global.discourse-cdn.com/julialang/original/3X/1/f/1f168106654f654695eac31a78e6e9b57881bd2b.png)

The code snippet I’ve included I checked up against the usual suspects here are the results for my machine:

```julia
@report_opt target_modules=(@ __MODULE__ ,)  
No errors detected

@benchmark PolyDataTemplate($save_location, $Points, $(["Kernel", "KernelGradient"]), $Kernel, $KernelGradient)
BenchmarkTools.Trial: 315 samples with 1 evaluation.
 Range (min … max): 13.444 ms … 21.689 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 15.613 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 15.856 ms ± 1.084 ms ┊ GC (mean ± σ): 0.28% ± 2.11%

               ▅▇ ▃▅█▇▁▄▃▁▂
  ▃▁▃▄▁▁▁▃▄▄▄▆▇██▇█████████▇▇▆█▇▆▄▆▄▃▃▃▄▄▃▃▃▄▃▁▁▃▁▁▃▁▁▃▁▃▁▁▁▃ ▄
  13.4 ms Histogram: frequency by time 19.9 ms <

 Memory estimate: 606.35 KiB, allocs estimate: 387.

```

I’ve only found through `@profview` some GC call I can’t seem to understand / remove:

![image](https://global.discourse-cdn.com/julialang/original/3X/8/4/8430bc41748660d2f86c212eb444df7cf0bdc223.png)

Any suggestions/help is appreciated 🙂

Kind regards

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [February 11, 2024, 2:39pm UTC](https://discourse.julialang.org/t/help-me-optimize-my-special-solution-to-save-vtp-vtk-file/110064/2 "2024-02-11T14:39:52Z")

</div>

(I have a comment, not really related to your question. I also plot/save similar data and found that these days I could make 3D particle plots much nicer/faster in Make with SSAO enabled (for proper feel of depth). That way one might even just use binary matrix formats for storing as paraview is not needed amymore. Anyway, don’t want to derail the topic!)

On the topic: What exactly is the issue you see? The amount of allocations might be reasonable given that it is an IO task…