diff --git a/sbmy_control/analysis/power_spectrum.py b/sbmy_control/analysis/power_spectrum.py index 7a92370..d97e24f 100644 --- a/sbmy_control/analysis/power_spectrum.py +++ b/sbmy_control/analysis/power_spectrum.py @@ -1,7 +1,7 @@ import numpy as np import os -kmin = 1e-1 +nfirst_kmodes = 10 kmax = 2e0 Nk = 50 AliasingCorr=False @@ -24,12 +24,17 @@ def crop_field(field, Ncrop): field.L2 = field.N2*d2 -def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None): +def get_power_spectrum(field, nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, G=None): from pysbmy.power import PowerSpectrum from pysbmy.fft import FourierGrid from pysbmy.correlations import get_autocorrelation if G is None: + default_k_modes = PowerSpectrum(field.L0,field.L1,field.L2,int(field.N0),int(field.N1),int(field.N2),).FourierGrid.k_modes + if kmax is None: + kmax = default_k_modes[-1] + if kmax > default_k_modes[-1]: + raise ValueError(f"kmax ({kmax}) is larger than the maximum k mode available in the field ({default_k_modes[-1]}).") G = FourierGrid( field.L0, field.L1, @@ -37,8 +42,8 @@ def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None): int(field.N0), int(field.N1), int(field.N2), - k_modes=np.concat([PowerSpectrum(field.L0,field.L1,field.L2,int(field.N0),int(field.N1),int(field.N2),).FourierGrid.k_modes[:10],np.logspace( - np.log10(kmin), + k_modes=np.concat([default_k_modes[:nfirst_kmodes],np.logspace( + np.log10(default_k_modes[nfirst_kmodes]), np.log10(kmax), Nk, )]), @@ -51,13 +56,18 @@ def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None): return G, k, Pk -def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk, G=None): +def get_cross_correlations(field_A, field_B, nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, G=None): from pysbmy.power import PowerSpectrum from pysbmy.fft import FourierGrid from pysbmy.correlations import get_crosscorrelation if G is None: + default_k_modes = PowerSpectrum(field_A.L0,field_A.L1,field_A.L2,int(field_A.N0),int(field_A.N1),int(field_A.N2),).FourierGrid.k_modes + if kmax is None: + kmax = default_k_modes[-1] + if kmax > default_k_modes[-1]: + raise ValueError(f"kmax ({kmax}) is larger than the maximum k mode available in the field ({default_k_modes[-1]}).") G = FourierGrid( field_A.L0, field_A.L1, @@ -65,8 +75,8 @@ def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk, G=None int(field_A.N0), int(field_A.N1), int(field_A.N2), - k_modes=np.concat([PowerSpectrum(field_A.L0,field_A.L1,field_A.L2,int(field_A.N0),int(field_A.N1),int(field_A.N2),).FourierGrid.k_modes[:10],np.logspace( - np.log10(kmin), + k_modes=np.concat([default_k_modes[:nfirst_kmodes],np.logspace( + np.log10(default_k_modes[nfirst_kmodes]), np.log10(kmax), Nk, )]), @@ -104,7 +114,7 @@ def plot_power_spectra(filenames, yticks = np.linspace(0.9,1.1,11), bound1=0.01, bound2=0.02, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, figsize=(8,4), @@ -139,7 +149,7 @@ def plot_power_spectra(filenames, color=colors[i], linestyle=linestyles[i], marker=markers[i],), - power_args=dict(kmin=kmin, + power_args=dict(nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk), ) @@ -178,7 +188,7 @@ def plot_cross_correlations(filenames_A, yticks = np.linspace(0.99,1.001,12), bound1=0.001, bound2=0.002, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, figsize=(8,4), @@ -217,7 +227,7 @@ def plot_cross_correlations(filenames_A, color=colors[i], linestyle=linestyles[i], marker=markers[i],), - power_args=dict(kmin=kmin, + power_args=dict(nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk), ) @@ -283,8 +293,15 @@ def console_main(): parser.add_argument('-yrp', '--ylim_power', type=float, nargs=2, default=[0.9,1.1], help='Y-axis limits.') parser.add_argument('-yrc', '--ylim_corr', type=float, nargs=2, default=[0.99,1.001], help='Y-axis limits.') parser.add_argument('--crop', type=int, default=None, help='Remove the outter N pixels of the fields.') + parser.add_argument('--nfirst_kmodes', type=float, default=10, help='First k modes from all available k modes to be used in the power spectrum before the geomspace part.') + parser.add_argument('--kmax', type=float, default=None, help='Maximum k value for the power spectrum geomspace part. If None, it will be set to the maximum k mode available in the field.') + parser.add_argument('--Nk', type=int, default=50, help='Number of k values for the power spectrum geomspace part.') args = parser.parse_args() + + nfirst_kmodes = args.nfirst_kmodes + kmax = args.kmax + Nk = args.Nk if not args.power_spectrum and not args.cross_correlation: print('You must choose between power_spectrum and cross_correlation.') @@ -298,7 +315,7 @@ def console_main(): from pysbmy.field import read_field F_ref = read_field(args.directory+args.reference) crop_field(F_ref, args.crop) - G, _, Pk_ref = get_power_spectrum(F_ref, kmin=kmin, kmax=kmax, Nk=Nk) + G, _, Pk_ref = get_power_spectrum(F_ref, nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk) else: Pk_ref = None G = None @@ -325,7 +342,7 @@ def console_main(): yticks = yticks_power, bound1=0.01, bound2=0.02, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, ax=axes[0], @@ -344,7 +361,7 @@ def console_main(): yticks = yticks_corr, bound1=0.001, bound2=0.002, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, ax=axes[1], @@ -371,7 +388,7 @@ def console_main(): yticks = yticks_power, bound1=0.01, bound2=0.02, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, Ncrop=args.crop, @@ -392,7 +409,7 @@ def console_main(): yticks = yticks_corr, bound1=0.001, bound2=0.002, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, Ncrop=args.crop, diff --git a/sbmy_control/low_level.py b/sbmy_control/low_level.py index 26054bf..b5382c0 100644 --- a/sbmy_control/low_level.py +++ b/sbmy_control/low_level.py @@ -64,10 +64,14 @@ def stdout_redirector(stream): try: # Create a temporary file and redirect stdout to it tfile = tempfile.TemporaryFile(mode="w+b") + _redirect_stdout(tfile.fileno()) + # Yield to caller, then redirect stdout back to the saved fd yield + _redirect_stdout(saved_stdout_fd) + # Copy contents of temporary file to the given stream tfile.flush() tfile.seek(0, io.SEEK_SET) @@ -215,9 +219,9 @@ def get_progress_from_logfile(filename): with open(filename, "r") as f: lines = f.readlines() for line in lines: - if " operation " in line: + if "Step " in line: try: - splitted_line = line.split(" operation ")[1] + splitted_line = line.split("Step ")[1] splitted_line = splitted_line.split(":")[0] current_operation = int(splitted_line.split("/")[0]) total_operations = int(splitted_line.split("/")[1]) @@ -237,7 +241,7 @@ def progress_bar_from_logfile(filename:str, desc:str="", verbose:int=1, **kwargs from tqdm import tqdm from time import sleep k=0 - limit=600 + limit=6000 update_interval=0.2 wait_until_file_exists(filename, verbose=verbose, limit=limit) current_operation, total_operations = get_progress_from_logfile(filename) diff --git a/sbmy_control/main.py b/sbmy_control/main.py index 8bb7ff0..0dd9afe 100644 --- a/sbmy_control/main.py +++ b/sbmy_control/main.py @@ -199,25 +199,26 @@ def main(parsed_args): def check_consistency(card_dict, mode): - ## Check consistency of EvolutionMode and ModulePMCOLA - if mode == "PM" or mode == "allPM": - if card_dict["EvolutionMode"] != 1: - raise ValueError(f"EvolutionMode is not 1: EvolutionMode={card_dict["EvolutionMode"]}") - if card_dict["ModulePMCOLA"] != 1: - raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}") - elif mode == "tCOLA" or mode == "alltCOLA": - if card_dict["EvolutionMode"] != 2: - raise ValueError(f"EvolutionMode is not 2: EvolutionMode={card_dict["EvolutionMode"]}") - if card_dict["ModulePMCOLA"] != 1: - raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}") - elif mode == "LPT": - if card_dict["ModulePMCOLA"] !=0: - raise ValueError(f"ModulePMCOLA is not 0: ModulePMCOLA={card_dict["ModulePMCOLA"]}") - elif mode == "pre_sCOLA" or mode == "post_sCOLA" or mode == "sCOLA" or mode=="allsCOLA": - if card_dict["EvolutionMode"] != 3: - raise ValueError(f"EvolutionMode is not 3: EvolutionMode={card_dict["EvolutionMode"]}") - if card_dict["ModulePMCOLA"] != 1: - raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}") + # ## Check consistency of EvolutionMode and ModulePMCOLA + # if mode == "PM" or mode == "allPM": + # if card_dict["EvolutionMode"] != 1: + # raise ValueError(f"EvolutionMode is not 1: EvolutionMode={card_dict["EvolutionMode"]}") + # if card_dict["ModulePMCOLA"] != 1: + # raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}") + # elif mode == "tCOLA" or mode == "alltCOLA": + # if card_dict["EvolutionMode"] != 2: + # raise ValueError(f"EvolutionMode is not 2: EvolutionMode={card_dict["EvolutionMode"]}") + # if card_dict["ModulePMCOLA"] != 1: + # raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}") + # elif mode == "LPT": + # if card_dict["ModulePMCOLA"] !=0: + # raise ValueError(f"ModulePMCOLA is not 0: ModulePMCOLA={card_dict["ModulePMCOLA"]}") + # elif mode == "pre_sCOLA" or mode == "post_sCOLA" or mode == "sCOLA" or mode=="allsCOLA": + # if card_dict["EvolutionMode"] != 3: + # raise ValueError(f"EvolutionMode is not 3: EvolutionMode={card_dict["EvolutionMode"]}") + # if card_dict["ModulePMCOLA"] != 1: + # raise ValueError(f"ModulePMCOLA is not 1: ModulePMCOLA={card_dict["ModulePMCOLA"]}") + pass def console_main(): diff --git a/sbmy_control/monofonic.py b/sbmy_control/monofonic.py index 134b55b..861f14b 100644 --- a/sbmy_control/monofonic.py +++ b/sbmy_control/monofonic.py @@ -32,7 +32,7 @@ def main_monofonic(parsed_args): print_message("Running monofonic in slurm mode.", 1, "monofonic", verbose=parsed_args.verbose) slurm_dict=parse_arguments_slurm(parsed_args) main_dict=parse_arguments_main(parsed_args) - slurm_script = slurm_dict["scripts"]+"monofonic.sh" + slurm_script = slurm_dict["scripts"]+f"monofonic_{main_dict['simname']}.sh" if not isfile(slurm_script) or parsed_args.force: print_message(f"SLURM script {slurm_script} does not exist (or forced). Creating it.", 2, "monofonic", verbose=parsed_args.verbose) diff --git a/sbmy_control/parameters_card.py b/sbmy_control/parameters_card.py index bbdc1a6..cfeda44 100644 --- a/sbmy_control/parameters_card.py +++ b/sbmy_control/parameters_card.py @@ -47,6 +47,15 @@ def register_arguments_card(parser:ArgumentParser): parser.add_argument("--OutputFCsDensity", type=str, default=None, help="Output FCs density file.") parser.add_argument("--OutputFCsSnapshot", type=str, default=None, help="Output FCs snapshot file.") parser.add_argument("--OutputRngStateLPT", type=str, default=None, help="Output RNG state file.") + ## Tests with phiBCs and density + parser.add_argument("--WriteGravPot", type=bool, default=False, help="Write gravitational potential.") + parser.add_argument("--OutputGravitationalPotentialBase", type=str, default=None, help="Output gravitational potential base.") + parser.add_argument("--MeshGravPot", type=int, default=None, help="Mesh for gravitational potential.") + parser.add_argument("--WriteDensity", type=bool, default=False, help="Write density.") + parser.add_argument("--OutputDensityBase", type=str, default=None, help="Output density base.") + parser.add_argument("--MeshDensity", type=int, default=None, help="Mesh for density.") + parser.add_argument("--LoadPhiBCs", type=bool, default=False, help="Load phiBCs.") + parser.add_argument("--InputPhiBCsBase", type=str, default=None, help="Input phiBCs file base.") parser.add_argument("--WriteReferenceFrame", type=bool, default=False, help="Write reference frame (COCA).") parser.add_argument("--OutputMomentaBase", type=str, default=None, help="Output momenta base (COCA).") parser.add_argument("--ReadReferenceFrame", type=bool, default=False, help="Read reference frame (COCA).") @@ -126,6 +135,16 @@ def parse_arguments_card(parsed_args): OutputFCsDensity=parsed_args.OutputFCsDensity, OutputFCsSnapshot=parsed_args.OutputFCsSnapshot, OutputRngStateLPT=parsed_args.OutputRngStateLPT, + ## Tests with phiBCs and density + WriteGravPot=parsed_args.WriteGravPot, + OutputGravitationalPotentialBase=parsed_args.OutputGravitationalPotentialBase, + MeshGravPot=parsed_args.MeshGravPot, + WriteDensity=parsed_args.WriteDensity, + OutputDensityBase=parsed_args.OutputDensityBase, + MeshDensity=parsed_args.MeshDensity, + LoadPhiBCs=parsed_args.LoadPhiBCs, + InputPhiBCsBase=parsed_args.InputPhiBCsBase, + ## Cosmological parameters h=cosmo_dict["h"], Omega_m=cosmo_dict["Omega_m"], Omega_b=cosmo_dict["Omega_b"], @@ -211,6 +230,17 @@ def parse_arguments_card(parsed_args): card_dict["OutputFCsSnapshot"] = main_dict["resultdir"]+"final_particles_"+main_dict["simname"]+".gadget3" if card_dict["OutputRngStateLPT"] is None: card_dict["OutputRngStateLPT"] = main_dict["workdir"]+"rng_state.h5" + ## Tests with phiBCs and density + if card_dict["OutputGravitationalPotentialBase"] is None: + card_dict["OutputGravitationalPotentialBase"] = main_dict["workdir"]+"gravpot_"+main_dict["simname"] + if card_dict["MeshGravPot"] is None: + card_dict["MeshGravPot"] = card_dict["N_PM_mesh"] + if card_dict["OutputDensityBase"] is None: + card_dict["OutputDensityBase"] = main_dict["workdir"]+"density_"+main_dict["simname"] + if card_dict["MeshDensity"] is None: + card_dict["MeshDensity"] = card_dict["N_PM_mesh"] + if card_dict["InputPhiBCsBase"] is None: + card_dict["InputPhiBCsBase"] = main_dict["workdir"]+"gravpot_tCOLA" if card_dict["OutputMomentaBase"] is None: card_dict["OutputMomentaBase"] = main_dict["workdir"]+"momenta_"+main_dict["simname"]+"_" if card_dict["InputMomentaBase"] is None: @@ -333,6 +363,15 @@ def create_parameter_card_dict( OutputFCsDensity:str = 'fcs_density.h5', OutputFCsSnapshot:str = 'fcs_particles.gadget3', + ## Tests with phiBCs and density + WriteGravPot:bool = True, + OutputGravitationalPotentialBase:str = 'gravitational_potential.h5', + MeshGravPot:int = 128, + WriteDensity:bool = False, + OutputDensityBase:str = 'density.h5', + MeshDensity:int = 128, + LoadPhiBCs:bool = False, + InputPhiBCsBase:str = 'gravitational_potential.h5', ## COCA parameters WriteReferenceFrame:bool = False, OutputMomentaBase:str = 'momenta_', @@ -403,6 +442,46 @@ def create_parameter_card_dict( OutputLPTPotential1=OutputLPTPotential1, OutputLPTPotential2=OutputLPTPotential2, OutputTilesBase=OutputTilesBase, + + # Tests with phiBCs and density + WriteGravPot=int(WriteGravPot), + OutputGravitationalPotentialBase=OutputGravitationalPotentialBase, + MeshGravPot=MeshGravPot, + WriteDensity=int(WriteDensity), + OutputDensityBase=OutputDensityBase, + MeshDensity=MeshDensity, + LoadPhiBCs=int(LoadPhiBCs), + InputPhiBCsBase=InputPhiBCsBase, + + + # P3M default parameters + n_Tiles = 60, + # Number of tiles in each dimension for the P3M code + nPairsForceDiagnostic =0, + # Number of particle-hole pairs to be used for the force diagnostic + nBinsForceDiagnostic =0, + maxTrialsForceDiagnostic =0, + OutputForceDiagnostic ='', + PrintOutputTimestepsLog = 0, + OutputTimestepsLog = "", + fac_dyn_custom = 0.03, + da_max_early_custom = 0.0013, + fac_H_custom = 0.05, + fac_bend = 0.3125, + sub_bend1 = 0.012, + sub_bend2 = 0.014, + fac_p3m_fit = 0.15, + da_max_late_custom = 0.02, + # Polynomial coefficients for the P3M time step limiter. Only used if use_p3m_fit=1 + use_p3m_fit = 1, + p3m_fit_coeff0 = -3.19544528, + p3m_fit_coeff1 = -0.01948223, + p3m_fit_coeff2 = 0.5180241, + p3m_fit_coeff3 = -1.08372045, + p3m_fit_coeff4 = -0.71808255, + p3m_fit_coeff5 = -0.14370312, + RunForceDiagnostic = 0, + h=h, Omega_m=Omega_m, Omega_b=Omega_b, diff --git a/sbmy_control/parameters_monofonic.py b/sbmy_control/parameters_monofonic.py index e4a2169..d9cf896 100644 --- a/sbmy_control/parameters_monofonic.py +++ b/sbmy_control/parameters_monofonic.py @@ -95,8 +95,6 @@ def get_config_from_dict(monofonic_dict): "Omega_L": monofonic_dict["Omega_q"], "H0": monofonic_dict["h"]*100., "n_s": monofonic_dict["n_s"], - "sigma_8": monofonic_dict["sigma8"], - "A_s": monofonic_dict["A_s"], "Tcmb": monofonic_dict["Tcmb"], "k_p": monofonic_dict["k_p"], "N_ur": monofonic_dict["N_ur"], @@ -110,7 +108,14 @@ def get_config_from_dict(monofonic_dict): "transfer": "CLASS", "ztarget": 2.5, "ZeroRadiation": False - } + } + + if "A_s" in monofonic_dict.keys() and monofonic_dict["A_s"] is not None: + config["cosmology"]["A_s"] = monofonic_dict["A_s"] + elif "sigma8" in monofonic_dict.keys() and monofonic_dict["sigma8"] is not None: + config["cosmology"]["sigma8"] = monofonic_dict["sigma8"] + else: + raise KeyError("Either A_s or sigma8 must be provided in the cosmology parameters.") config["random"] = { "generator": "NGENIC", diff --git a/sbmy_control/scripts/convert_snapshot_to_density.py b/sbmy_control/scripts/convert_snapshot_to_density.py index dadefed..f3935d9 100644 --- a/sbmy_control/scripts/convert_snapshot_to_density.py +++ b/sbmy_control/scripts/convert_snapshot_to_density.py @@ -1,9 +1,7 @@ -from pysbmy.density import get_density_pm_snapshot -from pysbmy.snapshot import read_snapshot import argparse -def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, 0.0, 0.0)): +def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, 0.0, 0.0), time=1.0): """ Convert a snapshot to a density field. @@ -18,6 +16,9 @@ def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, corner : tuple of float Corner of the box (x, y, z). """ + from pysbmy.density import get_density_pm_snapshot + from pysbmy.snapshot import read_snapshot + # Read the snapshot print("Reading snapshot...") snap = read_snapshot(snapshot_path) @@ -28,12 +29,86 @@ def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, # Calculate density print("Calculating density...") F = get_density_pm_snapshot(snap, N, N, N, corner[0], corner[1], corner[2]) + F.time = time # Set the time for the field # Write density to file print("Writing density...") F.write(output_path) print("Density written to", output_path) print("Done.") + + +def convert_snapshots_to_density(snapshot_base_path, output_path, N=None, corner=(0.0, 0.0, 0.0), time=1.0): + """ + Convert multiple snapshots to density fields. + + Parameters + ---------- + snapshot_base_path : str + Base path for the snapshot files. + output_path : str + Path to the output density file. + N : int + Size of the density field grid (N x N x N). + corner : tuple of float + Corner of the box (x, y, z). + """ + from pysbmy.density import get_density_pm, density_to_delta + from pysbmy.snapshot import read_snapshot + from pysbmy.field import Field + import numpy as np + + # Get all snapshot files with path "snapshot_base_path*" + import glob + + snapshot_files = glob.glob(snapshot_base_path + "*") + if not snapshot_files: + raise FileNotFoundError(f"No snapshot files found at {snapshot_base_path}") + + A = None + L0 = L1 = L2 = None + d0 = d1 = d2 = None + + for snapshot_path in snapshot_files: + # Read the snapshot + print(f"Reading snapshot {snapshot_path}...") + snap = read_snapshot(snapshot_path) + + if N is None: + N = snap.Np0 + + if L0 is None: + # Get the grid parameters from the snapshot + L0 = snap.L0 + L1 = snap.L1 + L2 = snap.L2 + d0 = L0 / N + d1 = L1 / N + d2 = L2 / N + else: + # Ensure the grid parameters match the first snapshot + if not (snap.L0 == L0 and snap.L1 == L1 and snap.L2 == L2): + raise ValueError(f"All snapshots must have the same grid parameters. Got {snap.L0}, {snap.L1}, {snap.L2} but expected {L0}, {L1}, {L2}.") + + if A is None: + # Initialize the density field + A = np.zeros((N, N, N), dtype=np.float32) + + # Calculate density + print("Calculating density for this snapshot...") + A = get_density_pm(snap.pos, A, d0=d0, d1=d1, d2=d2) + + # Convert density to delta field + print("Converting density to delta field...") + A = density_to_delta(A) + + # Create a Field object + F = Field(L0=L0,L1=L1,L2=L2,corner0=corner[0],corner1=corner[1],corner2=corner[2],rank=1,N0=N,N1=N,N2=N,time=time,data=A) + # Write density to file + print("Writing density...") + F.write(output_path) + print("Density written to", output_path) + print("Done.") def console_main(): @@ -70,13 +145,23 @@ def console_main(): args = parser.parse_args() - convert_snapshot_to_density( - snapshot_path=args.snapshot, - output_path=args.output, - N=args.N, - corner=args.corner, - ) + if args.snapshot.endswith("h5") or args.snapshot.endswith("hdf5") or args.snapshot.endswith("gadget3"): + # If the snapshot is a single file, convert it to density + convert_snapshot_to_density( + snapshot_path=args.snapshot, + output_path=args.output, + N=args.N, + corner=args.corner, + ) + else: + # If the snapshot is a base path, convert all snapshots to density + convert_snapshots_to_density( + snapshot_base_path=args.snapshot, + output_path=args.output, + N=args.N, + corner=args.corner, + ) if __name__ == "__main__": - console_main() \ No newline at end of file + console_main() diff --git a/sbmy_control/simbelmyne.py b/sbmy_control/simbelmyne.py index a106ae9..ac2a6ba 100644 --- a/sbmy_control/simbelmyne.py +++ b/sbmy_control/simbelmyne.py @@ -16,7 +16,7 @@ def main_simbelmyne(parsed_args): if isfile(log_file): # Remove the preexisting log file to allow for the progress_bar to be run normally from os import remove - oremove(log_file) + remove(log_file) command_args = ["simbelmyne", paramfile, log_file]