bugfix convert snapshots to density

sbmy_control/low_level.py

@@ -64,14 +64,10 @@ 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)
@@ -219,9 +215,9 @@ def get_progress_from_logfile(filename):
     with open(filename, "r") as f:
         lines = f.readlines()
     for line in lines:
-        if "Step " in line:
+        if " operation " in line:
             try:
-                splitted_line = line.split("Step ")[1]
+                splitted_line = line.split(" operation ")[1]
                 splitted_line = splitted_line.split(":")[0]
                 current_operation = int(splitted_line.split("/")[0])
                 total_operations = int(splitted_line.split("/")[1])
@@ -241,7 +237,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=6000
+    limit=600
     update_interval=0.2
     wait_until_file_exists(filename, verbose=verbose, limit=limit)
     current_operation, total_operations = get_progress_from_logfile(filename)

sbmy_control/main.py

@@ -199,26 +199,25 @@ 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"]}")
-    pass
+    ## 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"]}")
 
 
 def console_main():

sbmy_control/parameters_card.py

@@ -47,15 +47,6 @@ 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).")
@@ -135,16 +126,6 @@ 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"],
@@ -230,17 +211,6 @@ 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:
@@ -363,15 +333,6 @@ 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_',
@@ -442,46 +403,6 @@ 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,