Many bugs for long simulations and isolated simulations corrected

Author: Ezequiel92

src/analysis/compute_quantities/aggregators.jl

@@ -503,16 +503,17 @@ function integrateQty(
 
         if present_idx == 1
 
-            integrated_qty = 0.0u"Msun*yr^-1"
+            integrated_qty = 0.0u"Msun * yr^-1"
 
         else
 
             # Get the physical times
-            times = data_dict[:sim_data].snapshot_table[:, 5]
+            times = data_dict[:sim_data].snapshot_table[!, :physical_times]
+
             # Compute the time between snapshots
             Δt = times[present_idx] - times[present_idx - 1]
 
-            integrated_qty = sum(computeSFR(data_dict; age_resol=Δt); init=0.0u"Msun*yr^-1")
+            integrated_qty = sum(computeSFR(data_dict; age_resol=Δt); init=0.0u"Msun * yr^-1")
 
         end
 
@@ -531,7 +532,7 @@ function integrateQty(
         else
 
             # Get the physical times
-            times = data_dict[:sim_data].snapshot_table[:, 5]
+            times = data_dict[:sim_data].snapshot_table[!, :physical_times]
             # Compute the time between snapshots
             Δt = times[present_idx] - times[present_idx - 1]
 
@@ -666,19 +667,19 @@ function integrateQty(
 
     elseif quantity == :scale_factor
 
-        integrated_qty = data_dict[:sim_data].snapshot_table[data_dict[:snap_data].global_index, 3]
+        integrated_qty = data_dict[:sim_data].snapshot_table[data_dict[:snap_data].global_index, :scale_factors]
 
     elseif quantity == :redshift
 
-        integrated_qty = data_dict[:sim_data].snapshot_table[data_dict[:snap_data].global_index, 4]
+        integrated_qty = data_dict[:sim_data].snapshot_table[data_dict[:snap_data].global_index, :redshifts]
 
     elseif quantity == :physical_time
 
-        integrated_qty = data_dict[:sim_data].snapshot_table[data_dict[:snap_data].global_index, 5]
+        integrated_qty = data_dict[:sim_data].snapshot_table[data_dict[:snap_data].global_index, :physical_times]
 
     elseif quantity == :lookback_time
 
-        integrated_qty = data_dict[:sim_data].snapshot_table[data_dict[:snap_data].global_index, 6]
+        integrated_qty = data_dict[:sim_data].snapshot_table[data_dict[:snap_data].global_index, :lookback_times]
 
     elseif quantity == :ode_gas_it
 
@@ -1328,7 +1329,7 @@ function scatterQty(data_dict::Dict, quantity::Symbol)::Vector{<:Number}
         else
 
             # Get the physical times
-            times = data_dict[:sim_data].snapshot_table[:, 5]
+            times = data_dict[:sim_data].snapshot_table[!, :physical_times]
             # Compute the time between snapshots
             Δt = times[present_idx] - times[present_idx - 1]
 
@@ -1351,7 +1352,7 @@ function scatterQty(data_dict::Dict, quantity::Symbol)::Vector{<:Number}
         else
 
             # Get the physical times
-            times = data_dict[:sim_data].snapshot_table[:, 5]
+            times = data_dict[:sim_data].snapshot_table[!, :physical_times]
             # Compute the time between snapshots
             Δt = times[present_idx] - times[present_idx - 1]
 

src/analysis/compute_quantities/other.jl

@@ -114,7 +114,7 @@ function computeTimeTicks(
         scale_factors = [readTime(path) for path in snapshot_paths]
         redshifts = @. (1.0 / scale_factors) - 1.0
         physical_times = computeTime(scale_factors, readSnapHeader(first_snapshot))
-        lookback_times = last(physical_times) .- physical_times
+        lookback_times = maximum(physical_times) .- physical_times
 
     else
 
@@ -127,7 +127,7 @@ function computeTimeTicks(
         redshifts = zeros(length(snapshot_paths))
         # For non-cosmological simulations, the time in the snapshot is the physical time
         physical_times = [readTime(path) * u_time for path in snapshot_paths]
-        lookback_times = last(physical_times) .- physical_times
+        lookback_times = maximum(physical_times) .- physical_times
 
     end
 
@@ -260,7 +260,7 @@ function computeSFR(
     !isempty(ages) || return Unitful.MassFlow[]
 
     # Allocate memory
-    sfr = zeros(typeof(1.0u"Msun*yr^-1"), length(ages))
+    sfr = zeros(typeof(1.0u"Msun * yr^-1"), length(ages))
 
     # Find the stellar particles younger than `age_resol`
     idxs = map(x -> x <= age_resol, ages)

src/analysis/data_acquisition.jl

@@ -1257,11 +1257,13 @@ function makeSimulationTable(simulation_path::String)::DataFrame
         corresponding snapshot"))
     )
 
-    paths  = [snapshot_paths, groupcat_paths]
-    labels = [:snapshot_paths, :groupcat_paths]
-    rows   = [[1:length(snapshot_paths);], [1:length(groupcat_paths);]]
+    # Number of rows
+    n = length(snapshot_paths)
 
-    source_table = unstack(flatten(DataFrame(l=labels, p=paths, ids=rows), [:p, :ids]), :l, :p)
+    source_table = DataFrame(
+        snapshot_paths = snapshot_paths,
+        groupcat_paths = [groupcat_paths; fill(missing, n - length(groupcat_paths))],
+    )
 
     # Add the file name number column
     numbers = snap_source[:numbers]
@@ -1282,7 +1284,13 @@ function makeSimulationTable(simulation_path::String)::DataFrame
     # Add the lookback time column
     insertcols!(source_table, 6, :lookback_times => lookback_times; copycols=false)
 
-    return identity.(sort!(DataFrame(source_table), [:physical_times]))
+    # Sort the table by physical time
+    sort!(source_table, :physical_times)
+
+    # Add the row indices
+    insertcols!(source_table, 1, :row_id => 1:nrow(source_table))
+
+    return identity.(DataFrame(source_table))
 
 end
 
@@ -1938,9 +1946,9 @@ function internalUnits(quantity::String, path::String)::Union{Unitful.Quantity,U
 
         elseif dimensions == Unitful.𝐓
 
-            # From internal units to Myr, for non-cosmological simulations,
+            # From internal units to Myr for non-cosmological simulations,
             # and to a dimensionless quantity for cosmological simulations
-            return cosmological ? Unitful.NoUnits : IU.t_cosmo
+            return cosmological ? Unitful.NoUnits : IU.t_newton
 
         elseif dimensions == Unitful.𝐌 * Unitful.𝐋^-3
 

src/analysis/data_analysis.jl

@@ -4397,13 +4397,13 @@ function daEvolution(
 
     for (slice_index, sim_table_data) in pairs(iterator)
 
-        global_index  = sim_table_data[1]
-        scale_factor  = sim_table_data[3]
-        redshift      = sim_table_data[4]
-        physical_time = sim_table_data[5]
-        lookback_time = sim_table_data[6]
-        snapshot_path = sim_table_data[7]
-        groupcat_path = sim_table_data[8]
+        global_index  = sim_table_data[:row_id]
+        scale_factor  = sim_table_data[:scale_factors]
+        redshift      = sim_table_data[:redshifts]
+        physical_time = sim_table_data[:physical_times]
+        lookback_time = sim_table_data[:lookback_times]
+        snapshot_path = sim_table_data[:snapshot_paths]
+        groupcat_path = sim_table_data[:groupcat_paths]
 
         # Skip missing snapshots
         !ismissing(snapshot_path) || continue
@@ -4586,12 +4586,12 @@ function daVirialAccretion(
         :sim_data => sim_data,
         :snap_data => Snapshot(
             snapshot_path,
-            sim_table_data[1],
+            sim_table_data[:row_id],
             1,
-            sim_table_data[5],
-            sim_table_data[6],
-            sim_table_data[3],
-            sim_table_data[4],
+            sim_table_data[:physical_times],
+            sim_table_data[:lookback_times],
+            sim_table_data[:scale_factors],
+            sim_table_data[:redshifts],
             snapshot_header,
         ),
         :gc_data => GroupCatalog(groupcat_path, groupcat_header),
@@ -4620,13 +4620,13 @@ function daVirialAccretion(
 
     for (slice_index, sim_table_data) in pairs(iterator[2:end])
 
-        global_index  = sim_table_data[1]
-        scale_factor  = sim_table_data[3]
-        redshift      = sim_table_data[4]
-        physical_time = sim_table_data[5]
-        lookback_time = sim_table_data[6]
-        snapshot_path = sim_table_data[7]
-        groupcat_path = sim_table_data[8]
+        global_index  = sim_table_data[:row_id]
+        scale_factor  = sim_table_data[:scale_factors]
+        redshift      = sim_table_data[:redshifts]
+        physical_time = sim_table_data[:physical_times]
+        lookback_time = sim_table_data[:lookback_times]
+        snapshot_path = sim_table_data[:snapshot_paths]
+        groupcat_path = sim_table_data[:groupcat_paths]
 
         # Get the snapshot header
         snapshot_header = readSnapHeader(snapshot_path)
@@ -4810,12 +4810,12 @@ function daDiscAccretion(
         :sim_data => sim_data,
         :snap_data => Snapshot(
             snapshot_path,
-            sim_table_data[1],
+            sim_table_data[:row_id],
             1,
-            sim_table_data[5],
-            sim_table_data[6],
-            sim_table_data[3],
-            sim_table_data[4],
+            sim_table_data[:physical_times],
+            sim_table_data[:lookback_times],
+            sim_table_data[:scale_factors],
+            sim_table_data[:redshifts],
             snapshot_header,
         ),
         :gc_data => GroupCatalog(groupcat_path, groupcat_header),
@@ -4846,13 +4846,13 @@ function daDiscAccretion(
 
     for (slice_index, sim_table_data) in pairs(iterator[2:end])
 
-        global_index  = sim_table_data[1]
-        scale_factor  = sim_table_data[3]
-        redshift      = sim_table_data[4]
-        physical_time = sim_table_data[5]
-        lookback_time = sim_table_data[6]
-        snapshot_path = sim_table_data[7]
-        groupcat_path = sim_table_data[8]
+        global_index  = sim_table_data[:row_id]
+        scale_factor  = sim_table_data[:scale_factors]
+        redshift      = sim_table_data[:redshifts]
+        physical_time = sim_table_data[:physical_times]
+        lookback_time = sim_table_data[:lookback_times]
+        snapshot_path = sim_table_data[:snapshot_paths]
+        groupcat_path = sim_table_data[:groupcat_paths]
 
         # Get the snapshot header
         snapshot_header = readSnapHeader(snapshot_path)
@@ -4951,7 +4951,7 @@ function daSFRtxt(
     smooth::Int=0,
 )::NTuple{2,Vector{<:Number}}
 
-    snapshot_paths = filter(!ismissing, sim_data.snapshot_table[!, 7])
+    snapshot_paths = filter(!ismissing, sim_data.snapshot_table[!, :snapshot_paths])
 
     (
         !isempty(snapshot_paths) ||
@@ -5094,7 +5094,7 @@ function daCPUtxt(
     smooth::Int=0,
 )::NTuple{2,Vector{<:Number}}
 
-    snapshot_paths = filter(!ismissing, sim_data.snapshot_table[!, 7])
+    snapshot_paths = filter(!ismissing, sim_data.snapshot_table[!, :snapshot_paths])
 
     (
         !isempty(snapshot_paths) ||

src/analysis/filters.jl

@@ -1114,7 +1114,7 @@ function filterOldStars(data_dict::Dict)::Dict{Symbol,IndexType}
         end
 
         # Get the physical times
-        times = data_dict[:sim_data].snapshot_table[:, 5]
+        times = data_dict[:sim_data].snapshot_table[!, :physical_times]
 
         new_stars_idxs = map(t -> t > times[present_idx - 1], birth_times)
 

src/auxiliary_functions/histograms.jl

@@ -136,7 +136,7 @@ function histogram1D(
     width = (p_max - p_min) / n_bins
 
     # Allocate memory
-    histogram = zeros(eltype(values), n_bins)
+    histogram = zeros(typeof(first(values)), n_bins)
     counts = zeros(Int, n_bins)
 
     # Compute the histogram, ignoring NaNs and positions outside the grid range
@@ -227,7 +227,7 @@ function histogram1D(
     p_max = last(edges)
 
     # Allocate memory
-    histogram = zeros(eltype(values), n_bins)
+    histogram = zeros(typeof(first(values)), n_bins)
     counts = zeros(Int, n_bins)
 
     # Compute the histogram, ignoring NaNs and positions outside of range
@@ -420,7 +420,7 @@ function histogram2D(
     y_borders = (grid.y_bins[1] - h_bin_width, grid.y_bins[end] + h_bin_width)
 
     # Allocate memory
-    histogram = zeros(eltype(values), size(grid.grid))
+    histogram = zeros(typeof(first(values)), size(grid.grid))
     counts = zeros(Int, size(grid.grid))
 
     # Compute the histogram, ignoring NaNs and positions outside of range
@@ -529,7 +529,7 @@ function histogram2D(
     y_borders = (first(y_edges), last(y_edges))
 
     # Allocate memory
-    histogram = zeros(eltype(values), (n_x_bins, n_y_bins))
+    histogram = zeros(typeof(first(values)), (n_x_bins, n_y_bins))
     counts = zeros(Int, (n_x_bins, n_y_bins))
 
     # Compute the histogram, ignoring NaNs and positions outside of range
@@ -773,7 +773,7 @@ function histogram3D(
     z_borders = (grid.z_bins[1] - h_bin_width, grid.z_bins[end] + h_bin_width)
 
     # Allocate memory
-    histogram = zeros(eltype(values), size(grid.grid))
+    histogram = zeros(typeof(first(values)), size(grid.grid))
     counts = zeros(Int, size(grid.grid))
 
     # Compute the histogram, ignoring NaNs and positions outside of range

src/constants/globals.jl

@@ -660,6 +660,7 @@ Unit conversion factors.
   - `m_cosmo::Unitful.Mass`: From internal units of mass to ``\\mathrm{M_\\odot}``.
   - `t_cgs::Unitful.Time`: From internal units of time to ``\\mathrm{s}``.
   - `t_cosmo::Unitful.Time`: From internal units of time to ``\\mathrm{Myr}``.
+  - `t_newton::Unitful.Time`: From internal units (non-cosmological simulations) of time to ``\\mathrm{Myr}``.
   - `U_cgs::Unitful.Energy`: From internal units of specific energy to ``\\mathrm{erg \\, g^{-1}}``.
   - `rho_cgs::Unitful.Density`: From internal units of density to ``\\mathrm{g \\, cm^{-3}}``.
   - `P_Pa::Unitful.Pressure`: From internal units of pressure to ``\\mathrm{Pa}``.
@@ -678,6 +679,7 @@ struct InternalUnits
 
     t_cgs::Unitful.Time        # From internal units of time to s
     t_cosmo::Unitful.Time      # From internal units of time to Myr
+    t_newton::Unitful.Time     # From internal units (non-cosmological simulations) of time to Myr
 
     U_cgs::SpecificEnergy      # From internal units of specific energy to erg * g^-1
 
@@ -732,6 +734,7 @@ struct InternalUnits
         # Time conversion factors
         t_cgs = x_cgs / v_cgs
         t_cosmo = t_cgs |> u"Myr"
+        t_newton = 1.0u"Gyr" |> u"Myr"
 
         # Specific energy conversion factor
         U_cgs = v_unit^2 |> u"erg * g^-1"
@@ -753,6 +756,7 @@ struct InternalUnits
             m_cosmo,
             t_cgs,
             t_cosmo,
+            t_newton,
             U_cgs,
             rho_cgs,
             P_Pa,

src/plotting/convenience.jl

@@ -87,11 +87,14 @@ function snapshotReport(
             $(slice_n), the contents of $(simulation_path) are: \n$(simulation_table)"))
         )
 
+        # Select the ordinal index of the target snapshot
+        o_idx = findfirst(==(lpad(snap_number, 3, "0")), simulation_table.numbers)
+
         # Find the target row
-        snapshot_row = filter(:numbers => ==(lpad(snap_number, 3, "0")), simulation_table)
+        snapshot_row = simulation_table[o_idx, :]
 
         # Read the path to the target snapshot
-        snapshot_path = snapshot_row[1, :snapshot_paths]
+        snapshot_path = snapshot_row[:snapshot_paths]
 
         # Check that the snapshot path is not missing
         snapshot_filename = "\"$(SNAP_BASENAME)_$(lpad(snap_number, 3, "0"))\""
@@ -102,16 +105,13 @@ function snapshotReport(
         )
 
         # Read the path to the target group catalog file
-        groupcat_path = snapshot_row[1, :groupcat_paths]
+        groupcat_path = snapshot_row[:groupcat_paths]
 
         # Check if the simulation is cosmological
         cosmological = isSnapCosmological(snapshot_path)
 
         # Read the physical time since the Big Bang
-        physical_time = round(ustrip(u"Gyr", snapshot_row[1, :physical_times]), digits=2)
-
-        # Select the ordinal index of the target snapshot
-        o_idx = snapshot_row[1, :ids]
+        physical_time = round(ustrip(u"Gyr", snapshot_row[:physical_times]), digits=2)
 
         # Compute the number of snapshots in the folder
         snapshot_length = count(!ismissing, simulation_table[!, :snapshot_paths])
@@ -136,8 +136,8 @@ function snapshotReport(
 
         if cosmological
             # For cosmological simulations print the scale factor and the redshift
-            scale_factor = round(snapshot_row[1, :scale_factors], digits=3)
-            redshift = round(snapshot_row[1, :redshifts], digits=3)
+            scale_factor = round(snapshot_row[:scale_factors], digits=3)
+            redshift = round(snapshot_row[:redshifts], digits=3)
 
             println(file, "Scale factor:     $(scale_factor)")
             println(file, "Redshift:         $(redshift)")