minor changes

analysis/power_spectrum.py

@@ -15,10 +15,10 @@ def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None):
             field.L0,
             field.L1,
             field.L2,
-            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(
+            field.N0,
+            field.N1,
+            field.N2,
+            k_modes=np.concat([PowerSpectrum(field.L0,field.L1,field.L2,field.N0,field.N1,field.N2,).FourierGrid.k_modes[:10],np.logspace(
                 np.log10(kmin),
                 np.log10(kmax),
                 Nk,
@@ -43,10 +43,10 @@ def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk, G=None
             field_A.L0,
             field_A.L1,
             field_A.L2,
-            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(
+            field_A.N0,
+            field_A.N1,
+            field_A.N2,
+            k_modes=np.concat([PowerSpectrum(field_A.L0,field_A.L1,field_A.L2,field_A.N0,field_A.N1,field_A.N2,).FourierGrid.k_modes[:10],np.logspace(
                 np.log10(kmin),
                 np.log10(kmax),
                 Nk,

analysis/slices.py

@@ -25,7 +25,6 @@ def plot_imshow_with_reference( data_list,
     - cmap: colormap to be used for plotting
     """
     import matplotlib.pyplot as plt
-    from matplotlib import ticker
 
     if titles is None:
         titles = [None for f in data_list]
@@ -35,7 +34,7 @@ def plot_imshow_with_reference( data_list,
     elif isinstance(L, int) or isinstance(L, float):
         L = [L for data in data_list]
 
-    sep = 10 if L[0] < 50 else 20 if L[0] < 100 else 50 if L[0]<250 else 100 if L[0] < 500 else 200 if L[0] < 1000 else 500 if L[0] < 2500 else 1000
+    sep = 10 if L[0] < 50 else 20 if L[0] < 200 else 100
     ticks = [np.arange(0, l+1, sep)*len(dat)/l for l, dat in zip(L,data_list)]
     tick_labels = [np.arange(0, l+1, sep) for l in L]
 
@@ -66,8 +65,6 @@ def plot_imshow_with_reference( data_list,
             axes[0, i].set_xticklabels(tick_labels[i])
             axes[0, i].set_yticklabels(tick_labels[i])
             axes[0, i].set_xlabel('Mpc/h')
-            axes[0, i].xaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
-            axes[0, i].yaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
         fig.colorbar(im, ax=axes[0, :], orientation='vertical')
 
         # Plot the data compared to the reference
@@ -79,13 +76,11 @@ def plot_imshow_with_reference( data_list,
             axes[1, i].set_xticklabels(tick_labels[i])
             axes[1, i].set_yticklabels(tick_labels[i])
             axes[1, i].set_xlabel('Mpc/h')
-            axes[1, i].xaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
-            axes[1, i].yaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
         fig.colorbar(im, ax=axes[1, :], orientation='vertical')
         
         # Add the score on the plots
         for i, data in enumerate(data_list):
-            axes[1, i].text(0.5, 0.9, f"RMS: {score(data, reference):.2e}", fontsize=10, transform=axes[1, i].transAxes, color='white')
+            axes[1, i].text(0.5, 0.9, f"Score: {score(data, reference):.2e}", fontsize=10, transform=axes[1, i].transAxes, color='white')
         # plt.tight_layout()
     else:
 
@@ -98,8 +93,6 @@ def plot_imshow_with_reference( data_list,
             axes.set_xticklabels(tick_labels[0])
             axes.set_yticklabels(tick_labels[0])
             axes.set_xlabel('Mpc/h')
-            axes.xaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
-            axes.yaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
             fig.colorbar(im, ax=axes, orientation='vertical')
 
         else:
@@ -111,8 +104,6 @@ def plot_imshow_with_reference( data_list,
                 axes[i].set_xticklabels(tick_labels[i])
                 axes[i].set_yticklabels(tick_labels[i])
                 axes[i].set_xlabel('Mpc/h')
-                axes[i].xaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
-                axes[i].yaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
             fig.colorbar(im, ax=axes[:], orientation='vertical')
 
     return fig, axes
@@ -135,12 +126,12 @@ if __name__ == "__main__":
     parser.add_argument('-t', '--title', type=str, default=None, help='Title for the plot.')
     parser.add_argument('-log','--log_scale', action='store_true', help='Use log scale for the data.')
 
-    # register_arguments_cosmo(parser)
+    register_arguments_cosmo(parser)
 
     args = parser.parse_args()
 
     from pysbmy.field import read_field
-    # from pysbmy.cosmology import d_plus
+    from pysbmy.cosmology import d_plus
 
     ref_field = read_field(args.directory+args.reference) if args.reference is not None else None
     fields = [read_field(args.directory+f) for f in args.filenames]
@@ -180,6 +171,6 @@ if __name__ == "__main__":
     fig.suptitle(args.title)
     
     if args.output is not None:
-        fig.savefig(args.output,bbox_inches='tight')
+        fig.savefig(args.output)
     else:
-        fig.savefig(args.directory+'slices.jpg',bbox_inches='tight')
+        fig.savefig(args.directory+'slices.png')

convert_snapshot_to_density.py

@@ -1,54 +0,0 @@
-from pysbmy.density import get_density_pm_snapshot
-from pysbmy.snapshot import read_snapshot
-import argparse
-
-if __name__ == "__main__":
-    parser = argparse.ArgumentParser(description="Convert snapshot to density.")
-    parser.add_argument(
-        "-S",
-        "--snapshot",
-        type=str,
-        required=True,
-        help="Path to the snapshot file.",
-    )
-    parser.add_argument(
-        "-o",
-        "--output",
-        type=str,
-        required=True,
-        help="Path to the output density file.",
-    )
-    parser.add_argument(
-        "-N",
-        "--N",
-        type=int,
-        default=None,
-        help="Size of the density field grid (N x N x N).",
-    )
-    parser.add_argument(
-        "-c",
-        "--corner",
-        type=float,
-        nargs=3,
-        default=[0.0, 0.0, 0.0],
-        help="Corner of the box (x, y, z).",
-    )
-
-    args = parser.parse_args()
-
-    # Read the snapshot
-    print("Reading snapshot...")
-    snap = read_snapshot(args.snapshot)
-    
-    if args.N is None:
-        N = snap.Np0
-    else:
-        N = args.N
-
-    print("Calculating density...")
-    F=get_density_pm_snapshot(snap, N,N,N, args.corner[0],args.corner[1],args.corner[2])
-
-    print("Writing density...")
-    F.write(args.output)
-    print("Density written to", args.output)
-    print("Done.")

low_level.py

@@ -225,8 +225,6 @@ def get_progress_from_logfile(filename):
                 pass
         elif "Fatal" in line or "Error" in line:
             return -1, -1
-        elif "Everything done successfully, exiting." in line:
-            current_operation = total_operations
     return current_operation, total_operations
 
 
@@ -239,6 +237,7 @@ def progress_bar_from_logfile(filename:str, desc:str="", verbose:int=1, **kwargs
     k=0
     limit=600
     update_interval=0.2
+    sleep(2) # Wait for the process to be launched, and for the previous log file to be overwritten if necessary.
     wait_until_file_exists(filename, verbose=verbose, limit=limit)
     current_operation, total_operations = get_progress_from_logfile(filename)
     previous_operation = 0
@@ -267,4 +266,4 @@ def progress_bar_from_logfile(filename:str, desc:str="", verbose:int=1, **kwargs
     if k*update_interval >= limit:
         print_message(f"Progress bar timed out after {limit} seconds.", 3, "low level", verbose=verbose)
     
-    return
+    return

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.")
 
 
 def register_arguments_card_for_ICs(parser:ArgumentParser):
@@ -118,6 +127,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"],
@@ -203,6 +222,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"
     
     return card_dict
 
@@ -321,6 +351,16 @@ 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',
+
                     ## Cosmological parameters
                     h:float = 0.6732,
                     Omega_m:float = 0.302,
@@ -381,6 +421,17 @@ 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,
+
             h=h,
             Omega_m=Omega_m,
             Omega_b=Omega_b,

parameters_monofonic.py

@@ -71,7 +71,7 @@ def get_config_from_dict(monofonic_dict):
     config["setup"] = {
         "GridRes": monofonic_dict["gridres"],
         "BoxLength": monofonic_dict["boxlength"],
-        "zstart": 99.0,
+        "zstart": 999.0,
         "LPTorder": 2,
         "DoBaryons": False,
         "DoBaryonVrel": False,

scola.py

@@ -1,5 +1,3 @@
-
-
 def main_scola(parsed_args):
     from args_main import parse_arguments_main
     from low_level import print_starting_module, print_message, print_ending_module, progress_bar_from_logfile, wait_until_file_exists
@@ -52,16 +50,14 @@ def main_scola(parsed_args):
         print_message("sCOLA finished.", 1, "scola", verbose=parsed_args.verbose)
     
     elif parsed_args.execution == "slurm":
-        from slurm_submission import create_slurm_script, parse_arguments_slurm, limit_slurm_arrays
+        from slurm_submission import create_slurm_script, parse_arguments_slurm
         from args_main import parse_arguments_main
-        from parameters_card import parse_arguments_card
         import os
         print_message("Running scola in slurm mode.", 1, "scola", verbose=parsed_args.verbose)
         slurm_dict=parse_arguments_slurm(parsed_args)
         main_dict=parse_arguments_main(parsed_args)
-        card_dict=parse_arguments_card(parsed_args)
 
-        if (have_no_tiles(parsed_args) or parsed_args.force) and card_dict["N_tiles"]**3 < limit_slurm_arrays :
+        if have_no_tiles(parsed_args) or parsed_args.force:
             ## Submit all boxes
             print_message("Submitting all boxes.", 2, "scola", verbose=parsed_args.verbose)
             slurm_script = slurm_dict["scripts"]+"scola_"+main_dict["simname"]+".sh"
@@ -97,7 +93,6 @@ def main_scola(parsed_args):
 
         else:
             ## Submit missing boxes
-            from time import sleep
             missing_tiles_arrays = get_missing_tiles_arrays(parsed_args)
             print_message(f"Submitting missing boxes: {missing_tiles_arrays}.", 2, "scola", verbose=parsed_args.verbose)
 
@@ -132,7 +127,6 @@ def main_scola(parsed_args):
                     subprocess.run(command_args)
 
                 print_message("sCOLA job submitted.", 2, "scola", verbose=parsed_args.verbose)
-                sleep((missing_tiles[1]-missing_tiles[0])*1.0) # Sleep for a bit to avoid overloading the scheduler
 
                 os.remove(slurm_script) # Remove the script after submission (because it is specific to the missing tiles)
 
@@ -230,23 +224,16 @@ def have_no_tiles(parsed_args):
 def get_missing_tiles_arrays(parsed_args):
     from os.path import isfile
     from parameters_card import parse_arguments_card
-    from slurm_submission import limit_slurm_arrays
 
     card_dict = parse_arguments_card(parsed_args)
     nboxes_tot = int(parsed_args.N_tiles**3)
     missing_tiles_arrays = []
     in_sequence_of_missing = False
-    len_sequence_of_missing = 0
     for b in range(1,nboxes_tot+1):
         if not isfile(card_dict["OutputTilesBase"]+str(b)+".h5"):
-            len_sequence_of_missing += 1
             if not in_sequence_of_missing:
                 missing_tiles_arrays.append([b])
                 in_sequence_of_missing = True
-                len_sequence_of_missing = 1
-            if b%limit_slurm_arrays==limit_slurm_arrays-1:
-                missing_tiles_arrays[-1].append(b)
-                in_sequence_of_missing = False
         elif in_sequence_of_missing:
             missing_tiles_arrays[-1].append(b-1)
             in_sequence_of_missing = False
@@ -287,4 +274,4 @@ if __name__ == "__main__":
     if parsed_args.execution == "slurm":
         for b in range(1,nboxes_tot+1):
             wait_until_file_exists(f"{card_dict['OutputTilesBase']}{b}.h5", verbose=parsed_args.verbose, limit=5*60)
-    main_post_scola(parsed_args)
+    main_post_scola(parsed_args)

slurm_submission.py

@@ -2,7 +2,6 @@ from argparse import ArgumentParser
 from args_main import parse_arguments_main
 
 path_to_monofonic_binary = "/home/aubin/monofonic/build/monofonIC"
-limit_slurm_arrays=799
 
 def register_arguments_slurm(parser:ArgumentParser):
     """
@@ -134,16 +133,12 @@ def create_slurm_script(slurm_template:str,
          - simbelmyne
          - scola
      """
-    index_sub=0
     
     if array is not None and job != "scola":
         raise ValueError(f"Array job range provided for job type {job}.")
     if array is None and job == "scola":
         raise ValueError(f"Array job range not provided for job type {job}.")
-    elif job == "scola":
-        index_sub=array[0]//limit_slurm_arrays
-        array=(array[0]%limit_slurm_arrays, array[1]%limit_slurm_arrays)
-        
+    
     from os.path import isfile
     if not isfile(slurm_template):
         raise FileNotFoundError(f"SLURM template {slurm_template} does not exist.")
@@ -160,7 +155,7 @@ def create_slurm_script(slurm_template:str,
         case "simbelmyne":
             command_line = f"{job} {job_config_file} {job_log}"
         case "scola":
-            command_line = f"{job} {job_config_file} {job_log} "+f"-b $(({index_sub*limit_slurm_arrays} + $SLURM_ARRAY_TASK_ID))"
+            command_line = f"{job} {job_config_file} {job_log} "+"-b ${SLURM_ARRAY_TASK_ID}"
         case _:
             raise ValueError(f"Job type {job} not recognized.")
     
@@ -168,7 +163,7 @@ def create_slurm_script(slurm_template:str,
     with open(slurm_script, "w") as f:
         for line in template:
             if job_name is not None and "--job-name" in line:
-                line = f"#SBATCH --job-name={index_sub if index_sub!=0 else ''}{job_name}\n"
+                line = f"#SBATCH --job-name={job_name}\n"
             if array is not None and "--array" in line:
                 line = f"#SBATCH --array={array[0]}-{array[1]}\n"
             if array is not None and ("--output" in line  or "--error" in line):