bugfix convert snapshots to density

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)

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():

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,