Fft tools

Project.toml

@@ -8,20 +8,28 @@ ADIOS2 = "e0ce9d3b-0dbd-416f-8264-ccca772f60ec"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [compat]
 ADIOS2 = ">=1.4"
 CairoMakie = ">=0.13"
 DelimitedFiles = "1"
 Distributions = ">=0.25"
+FFTW = "1.9.0"
 Printf = "1"
 ProgressMeter = "1"
+Roots = "2"
 StatsBase = ">=0.34"
+TOML = "1.0.3"
 Unitful = "1"
+WriteVTK = "1.21.2"
 julia = ">=1.10"
 
 [extras]

src/EntityTools.jl

@@ -9,22 +9,38 @@ module EntityTools
     using ProgressMeter
     using ADIOS2
     using Printf
+    using TOML
+    using FFTW
+    using WriteVTK 
 
     include("IO/adios/adios.jl")
     include("IO/adios/fields.jl")
     include("IO/adios/particles.jl")
     include("IO/adios/spectra.jl")
+    include("IO/vtk/write_vtk.jl")
     include("calc/phase.jl")
     include("calc/spectra.jl")
+    include("calc/powerspectrum.jl")
     include("debug/parse_timing.jl")
+    include("units/structs.jl")
+    include("units/setup.jl")
+    include("units/length.jl")
+    include("units/speed.jl")
+    include("units/time.jl")
+    include("units/Bfield.jl")
+    include("units/temperature.jl")
 
     export phase_map,
             spectrum,
             parse_timing,
             EntityData,
+            EntityUnits,
             find_closest_time,
             read_particles,
             read_field,
-            read_spectrum
+            read_spectrum,
+            bp_to_vtk,
+            power_spectrum,
+            get_Temp,
 
 end

src/IO/adios/adios.jl

@@ -91,6 +91,39 @@ function read_times(sim_path, subfolder, read_steps=false)
     end
 end
 
+"""
+    construct_times(sim_path, subfolder, read_steps=false)
+
+Finds all available output times for a given `subfolder`.
+If `read_steps=true` is set it also returns an array of the output steps constructed from file names.
+"""
+function construct_times(sim_path, subfolder, read_steps=false)
+    # find all available files in the folder
+    files = readdir(joinpath(sim_path, subfolder))
+    Nfiles = length(files)
+    # read written steps
+    i_steps = Vector{Int64}(undef, Nfiles)
+    for i in 1:Nfiles 
+        i_steps[i] = parse(Int64, split(files[i], ".")[2])
+    end
+
+    # read run settings to get timestep
+    runsettings = read_run_settings(sim_path, subfolder)
+    dt = runsettings["algorithms.timestep.dt"]
+
+    # read times
+    times = Vector{Float64}(undef, Nfiles)
+    for i = 1:Nfiles
+        times[i] = i_steps[i] * dt
+    end
+
+    if read_steps
+        return times, i_steps
+    else
+        return times
+    end
+end
+
 """
     find_closest_time(run_info::EntityData, t::Real)
 
@@ -269,18 +302,28 @@ function read_run_settings(sim_path, subfolder)
 end
 
 """
-    EntityData(sim_path::String)
+    EntityData(sim_path::String; verbose=false)
 
 Reads an overview of the simulation into an `EntityData` struct.
 """
-function EntityData(sim_path::String)
+function EntityData(sim_path::String; verbose=false)
 
     if !isdir(sim_path)
         error("$sim_path not present!")
     end
 
     # check which output folders are present
-    subfolders = readdir(sim_path)
+    if verbose
+        @info "Checking for output folders in $sim_path"
+    end
+    all_files = readdir(sim_path)
+
+    # filter out only the directories
+    if verbose
+        @info "Filtering subfolders"
+    end
+    subfolders = filter(entry -> isdir(joinpath(sim_path, entry)), all_files)
+
     # make sure to not try to read a plots folder
     filter!(x -> x != "Plots", subfolders)
     filter!(x -> x != "plots", subfolders)
@@ -291,7 +334,10 @@ function EntityData(sim_path::String)
     end
 
     # reading output files and times
-    times, i_steps = read_times(sim_path, subfolders[1], true)
+    if verbose
+        @info "Reading output times and steps"
+    end
+    times, i_steps = construct_times(sim_path, subfolders[1], true)
 
     # if fields are written, read their info
     field_info = nothing
@@ -312,6 +358,9 @@ function EntityData(sim_path::String)
     end
 
     # read the settings of the simulation
+    if verbose
+        @info "Reading run settings"
+    end
     runsettings = read_run_settings(sim_path, subfolders[1])
 
     # construct an EntityData struct

src/IO/adios/fields.jl

@@ -18,18 +18,26 @@ function read_fields_reduced(run_info::EntityData, field_names::Vector{String};
         data[field_name] = Vector{Float64}(undef, length(i))
     end 
 
+    if verbose
+        P = Progress(length(i); desc="Reading and reducing fields:")
+    end
+
     # loop over files
     for (j, _step) in enumerate(i)
         # read the data of the current timestep
         step_data = read_field(run_info, field_names; 
                     i_step=run_info.steps[_step],
                     slice=slice,
-                    verbose=verbose)
+                    verbose=false)
 
         # reduce the fields
         for field_name in field_names
             data[field_name][j] = reduction_function(step_data[field_name])
         end
+
+        if verbose
+            next!(P)
+        end
     end
 
     return data

src/IO/adios/particles.jl

@@ -22,11 +22,13 @@ function read_particles_reduced(run_info::EntityData, species::Integer,
         data[prop] = Vector{Float64}(undef, length(i))
     end 
 
+    if verbose
+        P = Progress(length(i); desc="Reading and reducing particles:")
+    end
+
     # loop over files
     for (j, _step) in enumerate(i)
-        if verbose
-            @info "reading step $(run_info.steps[_step])..."
-        end
+
         # read the data of the current timestep
         step_data = read_particles(run_info, species, properties; 
                     i_step=run_info.steps[_step],
@@ -37,6 +39,10 @@ function read_particles_reduced(run_info::EntityData, species::Integer,
         for prop in properties
             data[prop][j] = reduction_function(step_data[prop])
         end    # remember to close the file
+
+        if verbose
+            next!(P)
+        end
     end
 
     return data

src/IO/vtk/write_vtk.jl

@@ -0,0 +1,25 @@
+"""
+    bp_to_vtk(data::EntityData, vtk_file::String, i_step::Integer)
+
+Writes the fields of a given time step from an EntityData object to a VTK file.
+"""
+function bp_to_vtk(data::EntityData, vtk_file::String, i_step::Integer)
+
+    dx = (data.settings["grid.extent"][2] - data.settings["grid.extent"][1]) / data.settings["grid.resolution"][1] + 1
+    x = data.settings["grid.extent"][1]:dx:data.settings["grid.extent"][2]
+
+    dy = (data.settings["grid.extent"][4] - data.settings["grid.extent"][3]) / data.settings["grid.resolution"][2] + 1
+    y = data.settings["grid.extent"][3]:dy:data.settings["grid.extent"][4]
+
+    dz = (data.settings["grid.extent"][6] - data.settings["grid.extent"][5]) / data.settings["grid.resolution"][3] + 1
+    z = data.settings["grid.extent"][5]:dz:data.settings["grid.extent"][6]
+
+    fields = read_field(data, data.fields.fields, i_step = i_step, verbose=true)
+
+    vtk_grid(vtk_file, x, y, z) do vtk
+        for field in data.fields.fields
+            println("writing $field")
+            vtk[field] = fields[field]
+        end
+    end
+end