diff --git a/plotter.py b/plotter.py index e3a1dc3..fadf93d 100644 --- a/plotter.py +++ b/plotter.py @@ -99,7 +99,7 @@ def quiver(ax, map_h, map_v, extent, key_v=None, lognorm=False, label="", **kwar hh, vv = gethv(map_h, map_v, extent) # get norm information - norm_v = np.sqrt(map_h ** 2 + map_v ** 2) + norm_v = np.sqrt(map_h**2 + map_v**2) max_v = np.max(norm_v) min_v = np.min(norm_v) @@ -456,7 +456,7 @@ class Plotter(Aggregator, BaseProcessor): run=run, **kwargs, ) - + if movie: filenames[run].append(plot_filename) @@ -968,8 +968,10 @@ class Plotter(Aggregator, BaseProcessor): if sinks: try: - self.current_processor.sinks() - data = pd.DataFrame(self.current_processor.get_value("/datasets/sinks")) + self.current_processor.load_sinks_rule() + data = pd.DataFrame( + self.current_processor.get_value("/datasets/load_sinks_rule") + ) part_pos = data[["x", "y", "z"]].values unit_length /= self.current_processor.lbox except KeyError: @@ -1435,7 +1437,7 @@ class Plotter(Aggregator, BaseProcessor): self.logger.warning( "Errorbar may be meaningless when ytransform is used" ) - + # Offset to start x when y in over a given value if wait_until_over is not None: offset = np.argmax(y > wait_until_over) @@ -1480,7 +1482,7 @@ class Plotter(Aggregator, BaseProcessor): (a, b, rho, _map_rule, stderr) = linregress(np.log10(x), np.log10(y)) self.logger.info( "Power law fit y = x^({}) * {} with R^2 = {} and error is {}".format( - a, 10 ** b, rho, stderr + a, 10**b, rho, stderr ) ) else: @@ -1504,12 +1506,12 @@ class Plotter(Aggregator, BaseProcessor): residual = errfunc(c, np.log10(x), np.log10(y), yerr / y) self.logger.info( "Power law fit y = x^({}) * {} with residual {}".format( - a, 10 ** b, residual + a, 10**b, residual ) ) if label is None: label = r"Power-law fit with index {:.1f}".format(a) - plt.plot(x, (10 ** b) * x ** a, label=label, **kwargs) + plt.plot(x, (10**b) * x**a, label=label, **kwargs) def _gen_from_log(self, logrule, name_y, name_x="time", description="Generated"): if name_x == "time": @@ -1723,7 +1725,7 @@ class Plotter(Aggregator, BaseProcessor): "/series/sinks_from_log/time", "/series/sinks_from_log/ssm", xunit=U.Myr, - yunit=U.Msun / U.pc ** 2, + yunit=U.Msun / U.pc**2, ), "Mass of the sinks as a function of time divided by surface", kind="run", diff --git a/snapshotprocessor.py b/snapshotprocessor.py index f022403..c143ce8 100644 --- a/snapshotprocessor.py +++ b/snapshotprocessor.py @@ -293,6 +293,37 @@ class SnapshotProcessor(HDF5Container): "id": U.none, "level": U.none, } + # Units sinks + unit_sinks = { + "id": U.none, + "msink": U.Msun, + "dmfsink": U.Msun, + "x": "unit_density", + "y": "unit_density", + "z": "unit_density", + "vx": "unit_velocity", + "vy": "unit_velocity", + "vz": "unit_velocity", + "rot_period": "unit_time", + "lx": {"unit_mass": 1, "unit_length": 2, "unit_time": -1}, + "ly": {"unit_mass": 1, "unit_length": 2, "unit_time": -1}, + "lz": {"unit_mass": 1, "unit_length": 2, "unit_time": -1}, + "acc_rate": {"unit_mass": 1, "unit_time": -1}, + "acc_lum": U.none, + "age": U.year, + "int_lum": U.none, + "Teff": U.K, + "level": U.none, + "mbh": "unit_mass", + "tform": "unit_time", + "del_mass": U.Msun, + "rho_gas": "unit_density", + "cs**2": {"unit_velocity": 2}, + "etherm": {"unit_mass": 1, "unit_velocity": 2}, + "vx_gas": "unit_velocity", + "vy_gas": "unit_velocity", + "vz_gas": "unit_velocity", + } G = 1.0 # Gravitational constant @@ -675,9 +706,10 @@ class SnapshotProcessor(HDF5Container): self.save_data(data, filename, group) return data - def save_data(self, data, filename, group, overwrite=False): + def save_data(self, data, filename, group, overwrite=False, units=None): self.logger.debug(f"Writing {filename}") hdf5 = tables.open_file(filename, mode="a") + units = units or self.unit_key try: nb_written = 0 for key in data: @@ -688,7 +720,7 @@ class SnapshotProcessor(HDF5Container): hdf5.create_array( f"/{group}", key, data[key], "", createparents=True ) - unit = self._get_units(self.unit_key[key]) + unit = self._get_units(units[key]) hdf5.get_node(f"/{group}/{key}").unit = unit nb_written += 1 else: @@ -775,10 +807,13 @@ class SnapshotProcessor(HDF5Container): def convert_hdf5(self, filename=None): self.load_destructured(save=False) + self.load_sinks() + sinks = {key: self.sinks[key].values for key in self.sinks} if filename is None: filename = self.snap_filename self.save_data(self.cells, filename, "cells") self.save_data(self.parts, filename, "parts") + self.save_data(sinks, filename, "sinks", units=self.unit_sinks) self.unload_destructured() def get_nml(self, nml_key): @@ -1600,12 +1635,15 @@ class SnapshotProcessor(HDF5Container): alpha_g = (2.0 / 3) * alpha_g return alpha_g - def _sinks(self): + def load_sinks(self): csv_name = f"{self.path}/output_{self.num:05}/sink_{self.num:05}.csv" if not os.path.exists(csv_name): # If ratarmount was used csv_name = f"{self.path}/output_{self.num:05}/output_{self.num:05}/sink_{self.num:05}.csv" + if not os.path.exists(csv_name): + self.sinks = {} + return {} f = open(csv_name) first_line = f.readline() second_line = f.readline() @@ -1642,6 +1680,7 @@ class SnapshotProcessor(HDF5Container): "Teff", ] df = pd.read_csv(csv_name, header=None, names=header) + self.sinks = df return {key: df[key].values for key in df} @@ -1935,39 +1974,10 @@ class SnapshotProcessor(HDF5Container): "/maps", dependencies=["coldens", "T_mwavg", "omega_mwavg"], ), - "sinks": Rule( - self._sinks, + "load_sinks_rule": Rule( + self.load_sinks, group="/datasets", - unit={ - "id": U.none, - "msink": U.Msun, - "dmfsink": U.Msun, - "x": "unit_density", - "y": "unit_density", - "z": "unit_density", - "vx": "unit_velocity", - "vy": "unit_velocity", - "vz": "unit_velocity", - "rot_period": "unit_time", - "lx": {"unit_mass": 1, "unit_length": 2, "unit_time": -1}, - "ly": {"unit_mass": 1, "unit_length": 2, "unit_time": -1}, - "lz": {"unit_mass": 1, "unit_length": 2, "unit_time": -1}, - "acc_rate": {"unit_mass": 1, "unit_time": -1}, - "acc_lum": U.none, - "age": U.year, - "int_lum": U.none, - "Teff": U.K, - "level": U.none, - "mbh": "unit_mass", - "tform": "unit_time", - "del_mass": U.Msun, - "rho_gas": "unit_density", - "cs**2": {"unit_velocity": 2}, - "etherm": {"unit_mass": 1, "unit_velocity": 2}, - "vx_gas": "unit_velocity", - "vy_gas": "unit_velocity", - "vz_gas": "unit_velocity", - }, + unit=self.unit_sinks, ), "pspec": Rule(self.pspec, "Power spectrum", "/hdf5"), "write_particles": Rule(self._write_particles, "Particles file", "/hdf5"), diff --git a/turbox/__init__.py b/turbox/__init__.py index ea45dc6..e25b1ff 100644 --- a/turbox/__init__.py +++ b/turbox/__init__.py @@ -2,6 +2,7 @@ from .turbox import ( get_pspec, get_pspec_slope, build_turbox_data, + rho_pdf, apply_rule_pdf, span_resolution, ) diff --git a/turbox/turbox_params.yml b/turbox/turbox_params.yml new file mode 100644 index 0000000..501c914 --- /dev/null +++ b/turbox/turbox_params.yml @@ -0,0 +1,85 @@ +plot : # Plot parameters + put_title : False # Add a title to plot + + # Maps + ntick : 6 # Number of ticks for maps + + # Overlays + vel_red : 40 # Take point each vel_red for velocities + + tight_layout : False # Wheter to apply tight_layout + time_fmt : "{:.3g} {}" # Time format string, 1st field is time and 2nd is unit + +rcParams: # Anything there will be passed to matplotlib's rcParams + {} # remove and replace by rcParams keys + + +pymses: # Parameters for Pymses reader + + # Source settings + variables : ["rho","vel"] # Read these grid variables + part_variables : [] + check_variables : False + order : '<' # Bit order + + + # Processing options + levelmax : 20 # Maximal AMR level visited when looking levels + fft : False # Quick and dirty rendering using FFT + verbose : False # Let pymses write on standart output + multiprocessing : True # Whether to use multiprocessing + + # Camera settings + center : [0.5, 0.5, 0.5] # Center of the image + zoom : 1. # Zoom of the image + map_size : 256 # Size of the computed maps in pixel + + # Filter parameters + filter : False # Enable filtering + min_coords : [0.35, 0.35, 0.45] + max_coords : [0.65, 0.65, 0.55] + +input: # Parameters on how to look for input files (= output from Ramses) + + log_prefix : "run.log" # Prefix of the log file + label_filename : "label.txt" # Name of the label file + nml_filename : "run.nml" # name of the namelist file + ramses_ism : False # If ramses-ism is used + +out: # Parameters for post processing + copy_info : True # Copy logs, nml and info files to the outdir + tag : "" # Tag for the image + interactive : True # Interactive mode (keep figures open) + save : True # Save the plots on the disk + ext : '.jpeg' # extension for plots + ext_subfolder : True # Separate production by extension in subfolders + fmt : "" # Format of the output filename for plots + # The following keys are accepted + # {out} : The output directory (where hdf5 files are also stored) + # {run} : Name of the relevant run + # {num} : Name of the input file (from Ramses) + # {ext} : Extension defined above (with the dot) + # {name} : Name of the rule + # {subfolder} : replaced by /{ext[1:]}/ if ext_subfolder is true + # {tag} : Tag defined above + # {nml[nml_key]} : The value of nml_key in the namelist (ex: {nml[amr_params/levelmin]}) + + +process: # General setting of the post-processor module + verbose : True # Give more infos on what is going on + num_process : 10 # Number of forks + save_cells : True # Save cells structure on disk + save_parts : True # Save particles on disk + unload_cells : True # Save memory usage + +rules: # Specific rules parameters + turb_energy_threshold : -1 # Remove invalid data (<0 = no threshold) + + +astrophysix: # Parameters for astrophysix and galactica + simu_fmt : "{tag}_{run}" # Format of the name of simulation + descr_fmt : "{tag}_{run}" # Format of the default description + # The following keys are accepted + # {run} : Name of the relevant run + # {tag} : Tag defined above + # {nml[nml_key]} : The value of nml_key in the namelist (ex: {nml[amr_params/levelmin]}) diff --git a/utils/runselector.py b/utils/runselector.py index d0c3c03..004377a 100644 --- a/utils/runselector.py +++ b/utils/runselector.py @@ -52,7 +52,6 @@ class RunSelector: time=None, unit_time=None, allow_nodata=False, - fallback_namelist=True, ): """ Select runs and outputs with several filter options. @@ -72,7 +71,7 @@ class RunSelector: "last" select only the last output. "all" preselect all outputs (default) - nml_filename : str name of the namelist (should be the same for all outputs) + nml_filename : str name of the default namelist (otherwise look into the output files) filter_name : str, filter runs by name. Default "*" filter_nml : tuple or list of tupple. @@ -107,15 +106,10 @@ class RunSelector: self.path_in = path_in self.nml_filename = nml_filename - self.fallback_nml = fallback_namelist self.allow_nodata = allow_nodata self.namelist = {} - do_tests = ( - (not self.fallback_nml) - or (len(filter_nml) > 0) - or (sort_run_by is not None) - ) + do_tests = not self.allow_nodata self.runs = self.get_runs( in_runs, filter_name, filter_nml, sort_run_by, do_tests=do_tests ) @@ -226,7 +220,16 @@ class RunSelector: def load_namelist(self, run, path=None): if path is None: - path = f"{self.path_in}/{run}/{self.nml_filename}" + names = glob.glob( + self.path_in + "/" + run + "/output_[0-9][0-9][0-9][0-9][0-9]" + ) + + i = 0 + path = self.path_in + "/" + run + "/" + self.nml_filename + while not os.path.exists(path) and i < len(names): + path = f"{names[i]}/namelist.txt" + i += 1 + return NamelistRecursive(f90nml.read(path)) def get_nml_value(self, nml_key, run): @@ -448,14 +451,6 @@ class RunSelector: else: return [] - # Be sure we have a namelist - if self.fallback_nml and run not in self.namelist: - self.logger.warning( - f"Used fallback namelist for run {run} from output {nums[0]}" - ) - path = f"{self.path_in}/{run}/output_{nums[0]:05}/namelist.txt" - self.namelist[run] = self.load_namelist(run, path=path) - # -- Time getter according to unit_time if unit_time is None: