From bd9b62e2bb0ec3db31f140e0e4566675e264f241 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Sat, 15 Nov 2025 10:04:34 -0600 Subject: [PATCH 01/12] construct units from output or parameter file --- Project.toml | 2 ++ src/EntityTools.jl | 10 +++++- src/units/setup.jl | 47 +++++++++++++++++++++++++++ src/units/structs.jl | 75 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 133 insertions(+), 1 deletion(-) create mode 100644 src/units/setup.jl create mode 100644 src/units/structs.jl diff --git a/Project.toml b/Project.toml index e0ff8f0..8288bc3 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] @@ -21,6 +22,7 @@ Distributions = ">=0.25" Printf = "1" ProgressMeter = "1" StatsBase = ">=0.34" +TOML = "1.0.3" Unitful = "1" julia = ">=1.10" diff --git a/src/EntityTools.jl b/src/EntityTools.jl index 01f4c3f..385eddd 100644 --- a/src/EntityTools.jl +++ b/src/EntityTools.jl @@ -9,6 +9,7 @@ module EntityTools using ProgressMeter using ADIOS2 using Printf + using TOML include("IO/adios/adios.jl") include("IO/adios/fields.jl") @@ -17,6 +18,12 @@ module EntityTools include("calc/phase.jl") include("calc/spectra.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") export phase_map, spectrum, @@ -24,6 +31,7 @@ module EntityTools EntityData, read_particles, read_field, - read_spectrum + read_spectrum, + EntityUnits end diff --git a/src/units/setup.jl b/src/units/setup.jl new file mode 100644 index 0000000..f8a6fe2 --- /dev/null +++ b/src/units/setup.jl @@ -0,0 +1,47 @@ +""" + EntityUnits(data::EntityData) + +Reads the normalized units from the first available file in the `fields`, `particles`, or `spectra` subfolders of the simulation path. +Returns an [`EntityUnits`](@ref) struct containing the normalized units for the simulation. +""" +function EntityUnits(data::EntityData) + + if !isnothing(data.fields) + file = select_first_file(data.path, "fields") + elseif !isnothing(data.particles) + file = select_first_file(data.path, "particles") + elseif !isnothing(data.spectra) + file = select_first_file(data.path, "spectra") + else + error("No file found in fields, particles, or spectra subfolders.") + end + + return EntityUnits( adios_attribute_data(file, "scales.skindepth0"), + adios_attribute_data(file, "scales.larmor0"), + adios_attribute_data(file, "scales.sigma0"), + adios_attribute_data(file, "scales.B0"), + adios_attribute_data(file, "scales.omegaB0"), + adios_attribute_data(file, "scales.q0"), + adios_attribute_data(file, "scales.V0"), + adios_attribute_data(file, "scales.n0") ) +end + +""" + EntityUnits(toml_path::String) + +Reads the normalized units from the simulations' parameter file at the specified path. +""" +function EntityUnits(toml_path::String) + if !isfile(toml_path) + error("File $toml_path not found!") + end + data = TOML.parsefile(toml_path) + extent = data["grid"]["extent"][1] + + return EntityUnits( data["scales"]["skindepth0"], + data["scales"]["larmor0"], + data["particles"]["ppc0"], + data["grid"]["extent"], + data["grid"]["resolution"], + data["grid"]["metric"]["metric"] ) +end \ No newline at end of file diff --git a/src/units/structs.jl b/src/units/structs.jl new file mode 100644 index 0000000..a830ff0 --- /dev/null +++ b/src/units/structs.jl @@ -0,0 +1,75 @@ +""" + EntityUnits{T <: Real} + +A struct to hold the normalized units for a plasma simulation. +The fields are: +- `skindepth0`: The skin depth in normalized units. +- `larmor0`: The Larmor radius in normalized units. +- `sigma0`: The magnetization parameter, calculated as (skindepth0/larmor0)^2. +- `B0`: The magnetic field strength in normalized units, calculated as 1/larmor0. +- `omegaB0`: The cyclotron frequency in normalized units, calculated as B0. +- `q0`: The charge normalization factor, calculated as 1/(n0 * skindepth0^2). +- `V0`: The volume normalization factor, calculated as dx0^size(extent, 1). +- `n0`: The number density normalization factor, calculated as ppc0 / V0. +""" +struct EntityUnits{T} + + skindepth0::T + larmor0::T + sigma0::T + B0::T + omegaB0::T + q0::T + V0::T + n0::T + + """ + EntityUnits(skindepth0::T, larmor0::T, sigma0::T, B0::T, omegaB0::T, q0::T, V0::T, n0::T) where T <: Real + + Constructs an `EntityUnits` struct from the provided parameters. + """ + EntityUnits(skindepth0::T, larmor0::T, sigma0::T, B0::T, omegaB0::T, q0::T, V0::T, n0::T) where T <: Real = new{T}(skindepth0, larmor0, sigma0, B0, omegaB0, q0, V0, n0) + + """ + EntityUnits(skindepth0::T, larmor0::T) where T <: Real + + Constructs an `EntityUnits` struct from the provided skin depth and Larmor radius. + The other fields are calculated based on these two parameters. + Charges and densities are set to 1. + """ + function EntityUnits(skindepth0::T, larmor0::T) where T <: Real + sigma0 = (skindepth0/larmor0)^2 + B0 = 1/larmor0 + new{T}(skindepth0, larmor0, sigma0, B0, B0, T(1), T(1), T(1)) + end + + """ + EntityUnits(skindepth0::T, larmor0::T, ppc0::T, extent::Vector{T}, resolution::Integer, metric::String) where T <: Real + + Constructs an `EntityUnits` struct from the provided parameters. + - `skindepth0`: The skin depth in normalized units. + - `larmor0`: The Larmor radius in normalized units. + - `ppc0`: The number of particles per cell. + - `extent`: The spatial extent of the simulation domain. + - `resolution`: The number of grid points in each dimension. + - `metric`: The metric of the simulation (currently only "minkowski" is supported). + + """ + function EntityUnits(skindepth0::T, larmor0::T, ppc0::T, extent::Vector{Any}, resolution::Vector{<:Integer}, metric::String) where T <: Real + + if metric != "minkowski" + error("Only minkowski metric is currently supported.") + end + + dx0 = (extent[1][2] - extent[1][1]) / resolution[1] + V0 = dx0^size(extent, 1) + n0 = ppc0 / V0 + q0 = 1/ (n0 * skindepth0^2) + sigma0 = (skindepth0/larmor0)^2 + B0 = 1/larmor0 + + new{T}(skindepth0, larmor0, sigma0, B0, B0, q0, V0, n0) + end + +end + From 920390a7fa5a0c45533e1db04255d98053382801 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Sat, 22 Nov 2025 22:04:37 -0600 Subject: [PATCH 02/12] initial commit --- src/units/Bfield.jl | 0 src/units/length.jl | 8 ++++++++ src/units/speed.jl | 8 ++++++++ src/units/time.jl | 7 +++++++ 4 files changed, 23 insertions(+) create mode 100644 src/units/Bfield.jl create mode 100644 src/units/length.jl create mode 100644 src/units/speed.jl create mode 100644 src/units/time.jl diff --git a/src/units/Bfield.jl b/src/units/Bfield.jl new file mode 100644 index 0000000..e69de29 diff --git a/src/units/length.jl b/src/units/length.jl new file mode 100644 index 0000000..08a9be5 --- /dev/null +++ b/src/units/length.jl @@ -0,0 +1,8 @@ +""" + rLi(unit::EntityUnits, mi, me) + +Larmor radius of a particle with mass mi in normalized units. +""" +function rLi(unit::EntityUnits, mi::T, me::T=T(1)) where T <: Real + return sqrt(2/unit.sigma0)*sqrt(mi/me)*unit.skindepth0 +end \ No newline at end of file diff --git a/src/units/speed.jl b/src/units/speed.jl new file mode 100644 index 0000000..9e4b254 --- /dev/null +++ b/src/units/speed.jl @@ -0,0 +1,8 @@ +""" + vA(unit::EntityUnits) + +Alfvén speed in normalized units. +""" +function vA(unit::EntityUnits) + return 1 / sqrt(unit.sigma0) +end \ No newline at end of file diff --git a/src/units/time.jl b/src/units/time.jl new file mode 100644 index 0000000..05c4dd0 --- /dev/null +++ b/src/units/time.jl @@ -0,0 +1,7 @@ +""" + +omega_p = sqrt(4*pi*n*q^2/m) +""" +function omega_i(data::EntityData) + +end \ No newline at end of file From 3e1638409e3cf68aa0c2663a08ff510d9cefb7fb Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Sun, 23 Nov 2025 04:17:42 +0000 Subject: [PATCH 03/12] simplification --- src/units/setup.jl | 28 +++++++++------------------- 1 file changed, 9 insertions(+), 19 deletions(-) diff --git a/src/units/setup.jl b/src/units/setup.jl index f8a6fe2..0fe7b19 100644 --- a/src/units/setup.jl +++ b/src/units/setup.jl @@ -1,29 +1,19 @@ """ EntityUnits(data::EntityData) -Reads the normalized units from the first available file in the `fields`, `particles`, or `spectra` subfolders of the simulation path. +Constructs the units from the simulation settings stored in `data`. Returns an [`EntityUnits`](@ref) struct containing the normalized units for the simulation. """ function EntityUnits(data::EntityData) - if !isnothing(data.fields) - file = select_first_file(data.path, "fields") - elseif !isnothing(data.particles) - file = select_first_file(data.path, "particles") - elseif !isnothing(data.spectra) - file = select_first_file(data.path, "spectra") - else - error("No file found in fields, particles, or spectra subfolders.") - end - - return EntityUnits( adios_attribute_data(file, "scales.skindepth0"), - adios_attribute_data(file, "scales.larmor0"), - adios_attribute_data(file, "scales.sigma0"), - adios_attribute_data(file, "scales.B0"), - adios_attribute_data(file, "scales.omegaB0"), - adios_attribute_data(file, "scales.q0"), - adios_attribute_data(file, "scales.V0"), - adios_attribute_data(file, "scales.n0") ) + return EntityUnits( data.settings["scales.skindepth0"], + data.settings["scales.larmor0"], + data.settings["scales.sigma0"], + data.settings["scales.B0"], + data.settings["scales.omegaB0"], + data.settings["scales.q0"], + data.settings["scales.V0"], + data.settings["scales.n0"] ) end """ From 8f213f84718785fa08e83d65aa179cc1ca3f095f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludwig=20B=C3=B6ss?= Date: Tue, 25 Nov 2025 15:35:21 +0100 Subject: [PATCH 04/12] first attempt --- Project.toml | 4 +- src/EntityTools.jl | 5 +- src/calc/powerspectrum.jl | 144 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 151 insertions(+), 2 deletions(-) create mode 100644 src/calc/powerspectrum.jl diff --git a/Project.toml b/Project.toml index e0ff8f0..f1c61d8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,13 +1,14 @@ name = "EntityTools" uuid = "b454753b-5441-4cda-8f2f-910cbcbb7d97" -version = "0.1.0" authors = ["Ludwig Böss and contributors"] +version = "0.1.0" [deps] 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" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" @@ -18,6 +19,7 @@ ADIOS2 = ">=1.4" CairoMakie = ">=0.13" DelimitedFiles = "1" Distributions = ">=0.25" +FFTW = "1.9.0" Printf = "1" ProgressMeter = "1" StatsBase = ">=0.34" diff --git a/src/EntityTools.jl b/src/EntityTools.jl index 659190b..1db191d 100644 --- a/src/EntityTools.jl +++ b/src/EntityTools.jl @@ -9,6 +9,7 @@ module EntityTools using ProgressMeter using ADIOS2 using Printf + using FFTW include("IO/adios/adios.jl") include("IO/adios/fields.jl") @@ -16,6 +17,7 @@ module EntityTools include("IO/adios/spectra.jl") include("calc/phase.jl") include("calc/spectra.jl") + include("calc/powerspectrum.jl") include("debug/parse_timing.jl") export phase_map, @@ -25,6 +27,7 @@ module EntityTools find_closest_time, read_particles, read_field, - read_spectrum + read_spectrum, + power_spectrum end diff --git a/src/calc/powerspectrum.jl b/src/calc/powerspectrum.jl new file mode 100644 index 0000000..918863d --- /dev/null +++ b/src/calc/powerspectrum.jl @@ -0,0 +1,144 @@ +""" + flatten_powerspec(dim::Int, data; verbose::Bool=false) + +Computes a 1D powerspectrum for a +""" +flatten_powerspec(dim::Int, data; verbose::Bool=false) = flatten_powerspec(Val(dim), data; verbose) + +""" + +""" +flatten_powerspec(::Val{1}, spectrum; verbose::Bool=false) = ifftshift(spectrum) + +""" + +""" +function flatten_powerspec(::Val{2}, spectrum; verbose::Bool=false) + N = size(spectrum, 1) + center = N ÷ 2 + 1 + + # Maximum possible radius is corner of the box, but usually we just go to Nyquist (N/2) + max_r = floor(Int, N / 2) + + # Arrays to store sum of magnitudes and count of pixels per radius bin + radial_sum = zeros(Float64, max_r + 1) + radial_count = zeros(Int, max_r + 1) + + if verbose + println("Computing radial average (1D spectrum)...") + p = Progress(N^2) + end + + @inbounds for j in 1:N, i in 1:N + # Calculate distance from center (DC component) + dx = i - center + dy = j - center + dist = sqrt(dx^2 + dy^2) + + # Determine integer bin index (0 to max_r) + bin_idx = round(Int, dist) + + # Only accumulate if within our interest range + if bin_idx <= max_r + # Julia arrays are 1-based, so map bin 0 -> index 1 + radial_sum[bin_idx + 1] += spectrum[i, j] + radial_count[bin_idx + 1] += 1 + end + + if verbose + # count up progress meter + next!(p) + end + end + + # avoid division by zero + radial_avg = zeros(Float64, length(radial_sum)) + @inbounds for i in 1:length(radial_sum) + if radial_count[i] > 0 + radial_avg[i] = radial_sum[i] / radial_count[i] + end + end + + return radial_avg +end + +""" + +""" +function flatten_powerspec(::Val{3}, spectrum; verbose::Bool=false) + N = size(spectrum, 1) + center = N ÷ 2 + 1 + + # Maximum possible radius is corner of the box, but usually we just go to Nyquist (N/2) + max_r = floor(Int, N / 2) + + # Arrays to store sum of magnitudes and count of pixels per radius bin + radial_sum = zeros(Float64, max_r + 1) + radial_count = zeros(Int, max_r + 1) + + if verbose + println("Computing radial average (1D spectrum)...") + p = Progress(N) + end + + @inbounds for k in 1:N, j in 1:N, i in 1:N + # Calculate distance from center (DC component) + dx = i - center + dy = j - center + dz = k - center + dist = sqrt(dx^2 + dy^2 + dz^2) + + # Determine integer bin index (0 to max_r) + bin_idx = round(Int, dist) + + # Only accumulate if within our interest range + if bin_idx <= max_r + # Julia arrays are 1-based, so map bin 0 -> index 1 + radial_sum[bin_idx + 1] += spectrum[i, j, k] + radial_count[bin_idx + 1] += 1 + end + + if verbose + # count up progress meter + next!(p) + end + end + + # avoid division by zero + radial_avg = zeros(Float64, length(radial_sum)) + for i in 1:length(radial_sum) + if radial_count[i] > 0 + radial_avg[i] = radial_sum[i] / radial_count[i] + end + end + + return radial_avg +end + +""" + +""" +function power_spectrum(field, box_width; verbose::Bool=false) + + if verbose + @info "Fourier transform of grid" + end + # fourier transform of grid + Fk = fftshift(fft(field)) + if verbose + @info "Fourier transform done" + @info "Computing powerspectrum" + end + # reduce from N-dim to 1-dim + Pk = flatten_powerspec(ndims(Fk), abs2.(Fk); verbose) + + if verbose + @info "Powerspectrum done" + end + # construct modes on grid + N = size(Fk, 1) + i_max = N ÷ 2 + k = [i/box_width for i = 0:i_max] + + return k, Pk +end \ No newline at end of file From eb58436d3ff9a6a9deda141e65ad8102d96d040e Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 3 Dec 2025 23:55:39 +0000 Subject: [PATCH 05/12] bugfix for new entity output structure --- src/IO/adios/adios.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/IO/adios/adios.jl b/src/IO/adios/adios.jl index ea70203..7f1f4f7 100644 --- a/src/IO/adios/adios.jl +++ b/src/IO/adios/adios.jl @@ -280,7 +280,11 @@ function EntityData(sim_path::String) end # check which output folders are present - subfolders = readdir(sim_path) + all_files = readdir(sim_path) + + # filter out only the directories + 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) From e2c8ab09c61b5a4d77c4983bf7963cbc229f4243 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 3 Dec 2025 23:56:10 +0000 Subject: [PATCH 06/12] added progress bar --- src/IO/adios/fields.jl | 10 +++++++++- src/IO/adios/particles.jl | 12 +++++++++--- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/src/IO/adios/fields.jl b/src/IO/adios/fields.jl index 9a771fb..6158378 100644 --- a/src/IO/adios/fields.jl +++ b/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 diff --git a/src/IO/adios/particles.jl b/src/IO/adios/particles.jl index d967bc0..25db3cc 100644 --- a/src/IO/adios/particles.jl +++ b/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 From 464c1f2a14dde6f1a7cd38bd7e28fbe03f78c184 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 10 Feb 2026 22:02:34 +0000 Subject: [PATCH 07/12] speed up IO --- src/IO/adios/adios.jl | 51 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 48 insertions(+), 3 deletions(-) diff --git a/src/IO/adios/adios.jl b/src/IO/adios/adios.jl index 7f1f4f7..98e415d 100644 --- a/src/IO/adios/adios.jl +++ b/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,20 +302,26 @@ 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 + 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 @@ -295,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 @@ -316,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 From 8fdf108f279d78c98af5eb8564945efab28770f2 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Tue, 24 Feb 2026 21:29:27 +0000 Subject: [PATCH 08/12] temperature computation based on Sridhar et al. (2021) --- src/units/temperature.jl | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 src/units/temperature.jl diff --git a/src/units/temperature.jl b/src/units/temperature.jl new file mode 100644 index 0000000..6f09955 --- /dev/null +++ b/src/units/temperature.jl @@ -0,0 +1,37 @@ +using Roots + +const A = 1.187 +const B = 1.251 +const C = 0.714 +const D = 0.936 + +γ_e(E_int) = (A + B*E_int) / ( C + D * E_int ) + +# Sridhar et al. (2021), Eq. C5 +E(E_int, γ_bar, Γ) = (γ_bar - Γ) * Γ / ( 1 + γ_e(E_int) * (Γ^2 - 1) ) - E_int + +""" + get_Temp(N::T, T00::T, V1::T, V2::T, V3::T) where T <: Real + +Compute the temperature based on Sridhar et al. (2021, https://arxiv.org/pdf/2107.00263), Appendix C. +It uses the number density `N`, total energy `T00`, and velocity components `V1`, `V2`, and `V3`. +The function uses a root-finding method to solve for the internal energy and then converts it to temperature using the relation T = (2/3) * E_int. +""" +function get_Temp(N::T, T00::T, V1::T, V2::T, V3::T) where T <: Real + + v_bulk = V1^2 + V2^2 + V3^2 + γ_bar = T00 / N + Γ = 1 / sqrt(1 - v_bulk) + + _E(E_int) = E(E_int, γ_bar, Γ) + Temp = 0.0 + try + Temp = find_zero(_E, (-1000000.0, 1000000.0), verbose=false) + catch + Temp = NaN + end + + # internal energy to temperature conversion + Temp *= 2/3 + return Temp +end \ No newline at end of file From 1a511f6bcfd3b385c758f85979eecd9cc2c3a3e1 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Sun, 5 Apr 2026 01:59:00 +0000 Subject: [PATCH 09/12] fix powerspec calculation --- src/calc/powerspectrum.jl | 333 ++++++++++++++++++++++++++------------ 1 file changed, 231 insertions(+), 102 deletions(-) diff --git a/src/calc/powerspectrum.jl b/src/calc/powerspectrum.jl index 918863d..aaf4fc0 100644 --- a/src/calc/powerspectrum.jl +++ b/src/calc/powerspectrum.jl @@ -1,144 +1,273 @@ +using FFTW + """ - flatten_powerspec(dim::Int, data; verbose::Bool=false) + power_spectrum(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; cells_per_skindepth = 1.0) -Computes a 1D powerspectrum for a +Computes the power spectrum of the given fields. +The fields are expected to be 3D arrays with dimensions (nx, ny, nz) or 2D arrays with dimensions (nx, ny). +The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. """ -flatten_powerspec(dim::Int, data; verbose::Bool=false) = flatten_powerspec(Val(dim), data; verbose) +function power_spectrum(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; cells_per_skindepth = 1.0) + if ndims(field_x) == 2 + return power_spectrum_2D(field_x, field_y, field_z; cells_per_skindepth) + elseif ndims(field_x) == 3 + return power_spectrum_3D(field_x, field_y, field_z; k_norm = cells_per_skindepth) + else + error("Unsupported number of dimensions: $(ndims(field_x)). Only 2D and 3D arrays are supported.") + end +end """ + power_spectrum(field::AbstractArray; cells_per_skindepth = 1.0) +Computes the power spectrum of a single scalar field. +The field is expected to be a 3D array with dimensions (nx, ny, nz) +or a 2D array with dimensions (nx, ny). +The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. """ -flatten_powerspec(::Val{1}, spectrum; verbose::Bool=false) = ifftshift(spectrum) +function power_spectrum(field::AbstractArray; cells_per_skindepth = 1.0) + if ndims(field) == 2 + return power_spectrum_2D(field; cells_per_skindepth = cells_per_skindepth) + elseif ndims(field) == 3 + return power_spectrum_3D(field; cells_per_skindepth = cells_per_skindepth) + else + error("Unsupported number of dimensions: $(ndims(field)). Only 2D and 3D arrays are supported.") + end +end + + +# ===================================================================== +# 2D Power Spectrum Implementation +# ===================================================================== """ + power_spectrum_2D(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; cells_per_skindepth = 1.0) +Computes the 2D power spectrum of the given fields. The fields are expected to be 3D arrays with dimensions (nx, ny, nz) or 2D arrays with dimensions (nx, ny). +The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. """ -function flatten_powerspec(::Val{2}, spectrum; verbose::Bool=false) - N = size(spectrum, 1) - center = N ÷ 2 + 1 +function power_spectrum_2D(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; cells_per_skindepth = 1.0) + nx, ny = size(field_x, 1), size(field_x, 2) + + # Calculate frequencies + kx = fftfreq(nx, cells_per_skindepth) + ky = rfftfreq(ny, cells_per_skindepth) + kmin = ky[2] - # Maximum possible radius is corner of the box, but usually we just go to Nyquist (N/2) - max_r = floor(Int, N / 2) + # Perform 2D real fourier transform + # (halves the second dimension as specified by the tuple (2,1)) + field_x_fft = rfft(field_x, (2, 1)) + field_y_fft = rfft(field_y, (2, 1)) + field_z_fft = rfft(field_z, (2, 1)) - # Arrays to store sum of magnitudes and count of pixels per radius bin - radial_sum = zeros(Float64, max_r + 1) - radial_count = zeros(Int, max_r + 1) + nx_fft, ny_fft = size(field_x_fft, 1), size(field_x_fft, 2) - if verbose - println("Computing radial average (1D spectrum)...") - p = Progress(N^2) - end + # Determine maximum k to setup bins correctly + kmax = sqrt(maximum(abs, kx)^2 + ky[end]^2) + k_bins = collect(0:kmin:(kmax - kmin/10)) + nbins = length(k_bins) - @inbounds for j in 1:N, i in 1:N - # Calculate distance from center (DC component) - dx = i - center - dy = j - center - dist = sqrt(dx^2 + dy^2) - - # Determine integer bin index (0 to max_r) - bin_idx = round(Int, dist) + amp = zeros(nbins) + + # Loop over Fourier modes and bin by kperp + for j in 1:ny_fft + # Multiply energy by 2 for internal RFFT frequencies + nf = (1 < j < ny_fft) ? 2.0 : 1.0 + ky_val = ky[j] - # Only accumulate if within our interest range - if bin_idx <= max_r - # Julia arrays are 1-based, so map bin 0 -> index 1 - radial_sum[bin_idx + 1] += spectrum[i, j] - radial_count[bin_idx + 1] += 1 - end - - if verbose - # count up progress meter - next!(p) + for i in 1:nx_fft + kx_val = kx[i] + kperp = sqrt(kx_val^2 + ky_val^2) + + # Compute bin index + bin_idx = floor(Int, kperp / kmin) + 1 + + if 1 <= bin_idx <= nbins + # Calculate energy on the fly + energy = (abs2(field_x_fft[i, j]) + + abs2(field_y_fft[i, j]) + + abs2(field_z_fft[i, j])) * nf + + amp[bin_idx] += energy + end end end + + return k_bins, amp +end + +""" + power_spectrum_2D(field::AbstractArray; cells_per_skindepth = 1.0) +Computes the 2D power spectrum of a single scalar field. The field is expected to be a 2D array with dimensions (nx, ny). The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. +""" +function power_spectrum_2D(field::AbstractArray; cells_per_skindepth = 1.0) + nx, ny = size(field) + + # Frequencies: rfft halves the dimension listed first in the region tuple. + # Below, rfft(field, (2, 1)) will halve dimension 2 (ny). + kx = fftfreq(nx, cells_per_skindepth) + ky = rfftfreq(ny, cells_per_skindepth) + kmin = ky[2] + + # 2D RFFT + field_fft = rfft(field, (2, 1)) + nx_fft, ny_fft = size(field_fft) - # avoid division by zero - radial_avg = zeros(Float64, length(radial_sum)) - @inbounds for i in 1:length(radial_sum) - if radial_count[i] > 0 - radial_avg[i] = radial_sum[i] / radial_count[i] + # Setup bins + kmax = sqrt(maximum(abs, kx)^2 + ky[end]^2) + k_bins = collect(0:kmin:(kmax - kmin/10)) + nbins = length(k_bins) + + amp = zeros(nbins) + + # Single-pass loop fusion (Column-major iteration) + @inbounds for j in 1:ny_fft + # Reality condition: double energy for internal frequencies + nf = (1 < j < ny_fft) ? 2.0 : 1.0 + ky_val = ky[j] + + for i in 1:nx_fft + kx_val = kx[i] + kperp = sqrt(kx_val^2 + ky_val^2) + + bin_idx = floor(Int, kperp / kmin) + 1 + + if 1 <= bin_idx <= nbins + energy = abs2(field_fft[i, j]) * nf + amp[bin_idx] += energy + end end end - return radial_avg + return k_bins, amp end +# ===================================================================== +# 3D Power Spectrum Implementation +# ===================================================================== + + """ + power_spectrum_3D(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; k_norm = 1.0) +Computes the 3D power spectrum of the given fields. The fields are expected to be 3D arrays with dimensions (nx, ny, nz). +The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. """ -function flatten_powerspec(::Val{3}, spectrum; verbose::Bool=false) - N = size(spectrum, 1) - center = N ÷ 2 + 1 +function power_spectrum_3D(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; cells_per_skindepth = 1.0) - # Maximum possible radius is corner of the box, but usually we just go to Nyquist (N/2) - max_r = floor(Int, N / 2) + nx, ny, nz = size(field_x) - # Arrays to store sum of magnitudes and count of pixels per radius bin - radial_sum = zeros(Float64, max_r + 1) - radial_count = zeros(Int, max_r + 1) + # Calculate frequencies for a 3D grid + kx = rfftfreq(nx, cells_per_skindepth) + ky = fftfreq(ny, cells_per_skindepth) + kz = fftfreq(nz, cells_per_skindepth) - if verbose - println("Computing radial average (1D spectrum)...") - p = Progress(N) - end + # Use the fundamental wavenumber for bin widths + kmin = kx[2] - @inbounds for k in 1:N, j in 1:N, i in 1:N - # Calculate distance from center (DC component) - dx = i - center - dy = j - center - dz = k - center - dist = sqrt(dx^2 + dy^2 + dz^2) - - # Determine integer bin index (0 to max_r) - bin_idx = round(Int, dist) - - # Only accumulate if within our interest range - if bin_idx <= max_r - # Julia arrays are 1-based, so map bin 0 -> index 1 - radial_sum[bin_idx + 1] += spectrum[i, j, k] - radial_count[bin_idx + 1] += 1 - end - - if verbose - # count up progress meter - next!(p) - end - end + # Perform 3D RFFT (defaults to transforming all 3 dimensions) + field_x_fft = rfft(field_x) + field_y_fft = rfft(field_y) + field_z_fft = rfft(field_z) - # avoid division by zero - radial_avg = zeros(Float64, length(radial_sum)) - for i in 1:length(radial_sum) - if radial_count[i] > 0 - radial_avg[i] = radial_sum[i] / radial_count[i] + nx_fft, ny_fft, nz_fft = size(field_x_fft) + + # Determine maximum k (the furthest corner of the 3D k-space box) + kmax = sqrt(kx[end]^2 + maximum(abs, ky)^2 + maximum(abs, kz)^2) + k_bins = collect(0:kmin:(kmax - kmin/10)) + nbins = length(k_bins) + + amp = zeros(nbins) + + # Precompute whether nx is even for the reality condition + nx_is_even = (nx % 2 == 0) + + # Single-pass loop fusion for 3D spherical shells + # Loop order corresponds to memory layout: innermost loop processes the 1st dimension. + @inbounds for k in 1:nz_fft + kz_val = kz[k] + + for j in 1:ny_fft + ky_val = ky[j] + + for i in 1:nx_fft + kx_val = kx[i] + + # 3D Wavenumber magnitude + k_mag = sqrt(kx_val^2 + ky_val^2 + kz_val^2) + + # Compute bin index mathematically + bin_idx = floor(Int, k_mag / kmin) + 1 + + if 1 <= bin_idx <= nbins + # Reality condition for rfft over the FIRST dimension: + # DC (i=1) and Nyquist (if nx is even and i=nx_fft) are unique. + # All other modes represent two conjugate frequencies, so we double their energy. + nf = (i == 1 || (nx_is_even && i == nx_fft)) ? 1.0 : 2.0 + + # Calculate energy on the fly + energy = (abs2(field_x_fft[i, j, k]) + + abs2(field_y_fft[i, j, k]) + + abs2(field_z_fft[i, j, k])) * nf + + amp[bin_idx] += energy + end + end end end - return radial_avg + return k_bins, amp end -""" """ -function power_spectrum(field, box_width; verbose::Bool=false) - - if verbose - @info "Fourier transform of grid" - end - # fourier transform of grid - Fk = fftshift(fft(field)) - if verbose - @info "Fourier transform done" - @info "Computing powerspectrum" - end - # reduce from N-dim to 1-dim - Pk = flatten_powerspec(ndims(Fk), abs2.(Fk); verbose) + power_spectrum_3D(field::AbstractArray; cells_per_skindepth = 1.0) +Computes the 3D power spectrum of a single scalar field. +The field is expected to be a 3D array with dimensions (nx, ny, nz). +The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. +""" +function power_spectrum_3D(field::AbstractArray; cells_per_skindepth = 1.0) + nx, ny, nz = size(field) + + # Frequencies: Default 3D rfft() halves the FIRST dimension (nx) + kx = rfftfreq(nx, cells_per_skindepth) + ky = fftfreq(ny, cells_per_skindepth) + kz = fftfreq(nz, cells_per_skindepth) + kmin = kx[2] + + # 3D RFFT + field_fft = rfft(field) + nx_fft, ny_fft, nz_fft = size(field_fft) - if verbose - @info "Powerspectrum done" + # Setup bins + kmax = sqrt(kx[end]^2 + maximum(abs, ky)^2 + maximum(abs, kz)^2) + k_bins = collect(0:kmin:(kmax - kmin/10)) + nbins = length(k_bins) + + amp = zeros(nbins) + nx_is_even = (nx % 2 == 0) + + # Single-pass loop fusion (Column-major iteration) + @inbounds for k in 1:nz_fft + kz_val = kz[k] + for j in 1:ny_fft + ky_val = ky[j] + for i in 1:nx_fft + kx_val = kx[i] + + k_mag = sqrt(kx_val^2 + ky_val^2 + kz_val^2) + bin_idx = floor(Int, k_mag / kmin) + 1 + + if 1 <= bin_idx <= nbins + # Reality condition for rfft over the FIRST dimension + nf = (i == 1 || (nx_is_even && i == nx_fft)) ? 1.0 : 2.0 + + energy = abs2(field_fft[i, j, k]) * nf + amp[bin_idx] += energy + end + end + end end - # construct modes on grid - N = size(Fk, 1) - i_max = N ÷ 2 - k = [i/box_width for i = 0:i_max] - - return k, Pk + + return k_bins, amp end \ No newline at end of file From 0a96c0c6ec03f7b9882de369e13891b8e871f443 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Sun, 5 Apr 2026 02:34:19 +0000 Subject: [PATCH 10/12] added option to convert a bp snapshot to avtk file --- Project.toml | 6 +++++- src/EntityTools.jl | 7 ++++++- src/IO/vtk/write_vtk.jl | 25 +++++++++++++++++++++++++ 3 files changed, 36 insertions(+), 2 deletions(-) create mode 100644 src/IO/vtk/write_vtk.jl diff --git a/Project.toml b/Project.toml index 6853104..aaeb2df 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "EntityTools" uuid = "b454753b-5441-4cda-8f2f-910cbcbb7d97" -authors = ["Ludwig Böss and contributors"] version = "0.1.0" +authors = ["Ludwig Böss and contributors"] [deps] ADIOS2 = "e0ce9d3b-0dbd-416f-8264-ccca772f60ec" @@ -11,9 +11,11 @@ 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" @@ -23,9 +25,11 @@ 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] diff --git a/src/EntityTools.jl b/src/EntityTools.jl index 1bb223d..c4dc7d4 100644 --- a/src/EntityTools.jl +++ b/src/EntityTools.jl @@ -11,11 +11,13 @@ module EntityTools 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") @@ -26,6 +28,7 @@ module EntityTools include("units/speed.jl") include("units/time.jl") include("units/Bfield.jl") + include("units/temperature.jl") export phase_map, spectrum, @@ -36,6 +39,8 @@ module EntityTools read_particles, read_field, read_spectrum, - power_spectrum + bp_to_vtk, + power_spectrum, + get_Temp, end diff --git a/src/IO/vtk/write_vtk.jl b/src/IO/vtk/write_vtk.jl new file mode 100644 index 0000000..e3fd6f9 --- /dev/null +++ b/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 From 9c9e440874e3a9b10810a9997fda4cbb1e046d4e Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 8 Apr 2026 04:13:40 +0000 Subject: [PATCH 11/12] minor updates --- src/calc/powerspectrum.jl | 55 +++++++++++++++++++++++++++++++-------- 1 file changed, 44 insertions(+), 11 deletions(-) diff --git a/src/calc/powerspectrum.jl b/src/calc/powerspectrum.jl index aaf4fc0..5764e6d 100644 --- a/src/calc/powerspectrum.jl +++ b/src/calc/powerspectrum.jl @@ -7,11 +7,12 @@ Computes the power spectrum of the given fields. The fields are expected to be 3D arrays with dimensions (nx, ny, nz) or 2D arrays with dimensions (nx, ny). The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. """ -function power_spectrum(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; cells_per_skindepth = 1.0) +function power_spectrum(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; + cells_per_skindepth = 1.0, verbose::Bool = false) if ndims(field_x) == 2 - return power_spectrum_2D(field_x, field_y, field_z; cells_per_skindepth) + return power_spectrum_2D(field_x, field_y, field_z; cells_per_skindepth, verbose) elseif ndims(field_x) == 3 - return power_spectrum_3D(field_x, field_y, field_z; k_norm = cells_per_skindepth) + return power_spectrum_3D(field_x, field_y, field_z; cells_per_skindepth, verbose) else error("Unsupported number of dimensions: $(ndims(field_x)). Only 2D and 3D arrays are supported.") end @@ -25,11 +26,11 @@ The field is expected to be a 3D array with dimensions (nx, ny, nz) or a 2D array with dimensions (nx, ny). The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. """ -function power_spectrum(field::AbstractArray; cells_per_skindepth = 1.0) +function power_spectrum(field::AbstractArray; cells_per_skindepth = 1.0, verbose::Bool = false) if ndims(field) == 2 - return power_spectrum_2D(field; cells_per_skindepth = cells_per_skindepth) + return power_spectrum_2D(field; cells_per_skindepth, verbose) elseif ndims(field) == 3 - return power_spectrum_3D(field; cells_per_skindepth = cells_per_skindepth) + return power_spectrum_3D(field; cells_per_skindepth, verbose) else error("Unsupported number of dimensions: $(ndims(field)). Only 2D and 3D arrays are supported.") end @@ -46,14 +47,20 @@ end Computes the 2D power spectrum of the given fields. The fields are expected to be 3D arrays with dimensions (nx, ny, nz) or 2D arrays with dimensions (nx, ny). The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. """ -function power_spectrum_2D(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; cells_per_skindepth = 1.0) +function power_spectrum_2D(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; + cells_per_skindepth::Real = 1.0, verbose::Bool = false) nx, ny = size(field_x, 1), size(field_x, 2) + verbose && @info "Calculating 2D power spectrum with dimensions: ($nx, $ny)" + + verbose && @info "Calculating frequencies for 2D grid..." # Calculate frequencies kx = fftfreq(nx, cells_per_skindepth) ky = rfftfreq(ny, cells_per_skindepth) kmin = ky[2] + + verbose && @info "Performing 2D RFFT on fields..." # Perform 2D real fourier transform # (halves the second dimension as specified by the tuple (2,1)) field_x_fft = rfft(field_x, (2, 1)) @@ -68,6 +75,8 @@ function power_spectrum_2D(field_x::AbstractArray, field_y::AbstractArray, field nbins = length(k_bins) amp = zeros(nbins) + + P = Progress(ny_fft, desc="Binning Fourier modes...") # Loop over Fourier modes and bin by kperp for j in 1:ny_fft @@ -91,6 +100,8 @@ function power_spectrum_2D(field_x::AbstractArray, field_y::AbstractArray, field amp[bin_idx] += energy end end + + verbose && next!(P) end return k_bins, amp @@ -100,15 +111,19 @@ end power_spectrum_2D(field::AbstractArray; cells_per_skindepth = 1.0) Computes the 2D power spectrum of a single scalar field. The field is expected to be a 2D array with dimensions (nx, ny). The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. """ -function power_spectrum_2D(field::AbstractArray; cells_per_skindepth = 1.0) +function power_spectrum_2D(field::AbstractArray; cells_per_skindepth::Real = 1.0, verbose::Bool = false) nx, ny = size(field) + verbose && @info "Calculating 2D power spectrum with dimensions: ($nx, $ny)" + + verbose && @info "Calculating frequencies for 2D grid..." # Frequencies: rfft halves the dimension listed first in the region tuple. # Below, rfft(field, (2, 1)) will halve dimension 2 (ny). kx = fftfreq(nx, cells_per_skindepth) ky = rfftfreq(ny, cells_per_skindepth) kmin = ky[2] + verbose && @info "Performing 2D RFFT on field..." # 2D RFFT field_fft = rfft(field, (2, 1)) nx_fft, ny_fft = size(field_fft) @@ -119,7 +134,9 @@ function power_spectrum_2D(field::AbstractArray; cells_per_skindepth = 1.0) nbins = length(k_bins) amp = zeros(nbins) - + + P = Progress(ny_fft, desc="Binning Fourier modes...") + # Single-pass loop fusion (Column-major iteration) @inbounds for j in 1:ny_fft # Reality condition: double energy for internal frequencies @@ -137,6 +154,8 @@ function power_spectrum_2D(field::AbstractArray; cells_per_skindepth = 1.0) amp[bin_idx] += energy end end + + verbose && next!(P) end return k_bins, amp @@ -153,10 +172,13 @@ end Computes the 3D power spectrum of the given fields. The fields are expected to be 3D arrays with dimensions (nx, ny, nz). The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. """ -function power_spectrum_3D(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; cells_per_skindepth = 1.0) +function power_spectrum_3D(field_x::AbstractArray, field_y::AbstractArray, field_z::AbstractArray; + cells_per_skindepth::Real = 1.0, verbose::Bool = false) nx, ny, nz = size(field_x) + verbose && @info "Calculating 3D power spectrum with dimensions: ($nx, $ny, $nz)" + verbose && @info "Calculating frequencies for 3D grid..." # Calculate frequencies for a 3D grid kx = rfftfreq(nx, cells_per_skindepth) ky = fftfreq(ny, cells_per_skindepth) @@ -165,6 +187,7 @@ function power_spectrum_3D(field_x::AbstractArray, field_y::AbstractArray, field # Use the fundamental wavenumber for bin widths kmin = kx[2] + verbose && @info "Performing 3D RFFT on fields..." # Perform 3D RFFT (defaults to transforming all 3 dimensions) field_x_fft = rfft(field_x) field_y_fft = rfft(field_y) @@ -181,6 +204,8 @@ function power_spectrum_3D(field_x::AbstractArray, field_y::AbstractArray, field # Precompute whether nx is even for the reality condition nx_is_even = (nx % 2 == 0) + + P = Progress(nz_fft, desc="Binning Fourier modes...") # Single-pass loop fusion for 3D spherical shells # Loop order corresponds to memory layout: innermost loop processes the 1st dimension. @@ -214,6 +239,8 @@ function power_spectrum_3D(field_x::AbstractArray, field_y::AbstractArray, field end end end + + verbose && next!(P) end return k_bins, amp @@ -226,15 +253,19 @@ Computes the 3D power spectrum of a single scalar field. The field is expected to be a 3D array with dimensions (nx, ny, nz). The `cells_per_skindepth` parameter is used to scale the frequencies appropriately. """ -function power_spectrum_3D(field::AbstractArray; cells_per_skindepth = 1.0) +function power_spectrum_3D(field::AbstractArray; cells_per_skindepth::Real = 1.0, verbose::Bool = false) nx, ny, nz = size(field) + + verbose && @info "Calculating 3D power spectrum with dimensions: ($nx, $ny, $nz)" + verbose && @info "Calculating frequencies for 3D grid..." # Frequencies: Default 3D rfft() halves the FIRST dimension (nx) kx = rfftfreq(nx, cells_per_skindepth) ky = fftfreq(ny, cells_per_skindepth) kz = fftfreq(nz, cells_per_skindepth) kmin = kx[2] + verbose && @info "Performing 3D RFFT on field..." # 3D RFFT field_fft = rfft(field) nx_fft, ny_fft, nz_fft = size(field_fft) @@ -247,6 +278,7 @@ function power_spectrum_3D(field::AbstractArray; cells_per_skindepth = 1.0) amp = zeros(nbins) nx_is_even = (nx % 2 == 0) + P = Progress(nz_fft, desc="Binning Fourier modes...") # Single-pass loop fusion (Column-major iteration) @inbounds for k in 1:nz_fft kz_val = kz[k] @@ -267,6 +299,7 @@ function power_spectrum_3D(field::AbstractArray; cells_per_skindepth = 1.0) end end end + verbose && next!(P) end return k_bins, amp From 21e9f4aaed920feab0b7403647e668f16e60f6e7 Mon Sep 17 00:00:00 2001 From: LudwigBoess Date: Wed, 8 Apr 2026 04:14:00 +0000 Subject: [PATCH 12/12] minor fix --- src/units/temperature.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/units/temperature.jl b/src/units/temperature.jl index 6f09955..0a969a7 100644 --- a/src/units/temperature.jl +++ b/src/units/temperature.jl @@ -11,13 +11,13 @@ const D = 0.936 E(E_int, γ_bar, Γ) = (γ_bar - Γ) * Γ / ( 1 + γ_e(E_int) * (Γ^2 - 1) ) - E_int """ - get_Temp(N::T, T00::T, V1::T, V2::T, V3::T) where T <: Real + get_Temp(N::T, T00::T, V1::T, V2::T, V3::T; adiabatic_index::Real=5/3) where T <: Real Compute the temperature based on Sridhar et al. (2021, https://arxiv.org/pdf/2107.00263), Appendix C. It uses the number density `N`, total energy `T00`, and velocity components `V1`, `V2`, and `V3`. The function uses a root-finding method to solve for the internal energy and then converts it to temperature using the relation T = (2/3) * E_int. """ -function get_Temp(N::T, T00::T, V1::T, V2::T, V3::T) where T <: Real +function get_Temp(N::T, T00::T, V1::T, V2::T, V3::T; adiabatic_index::Real=5/3) where T <: Real v_bulk = V1^2 + V2^2 + V3^2 γ_bar = T00 / N @@ -32,6 +32,6 @@ function get_Temp(N::T, T00::T, V1::T, V2::T, V3::T) where T <: Real end # internal energy to temperature conversion - Temp *= 2/3 + Temp *= (adiabatic_index - 1) return Temp end \ No newline at end of file