SPH for CSiBORG1 support (#109)

* Add edited files

* Update scripts

* Updates

* Little fix

* Little ifx

scripts_independent/field_sph_gadget.py

@@ -16,11 +16,11 @@
 Script to construct the density and velocity fields for a simulation snapshot.
 The SPH filter is implemented in the cosmotool package.
 """
-from argparse import ArgumentParser
-from os import remove
-from os.path import join, exists
 import subprocess
+from argparse import ArgumentParser
 from datetime import datetime
+from os import remove
+from os.path import exists, join
 
 import hdf5plugin  # noqa
 import numpy as np
@@ -137,7 +137,7 @@ def main(snapshot_path, output_path, resolution, scratch_space, SPH_executable,
     SPH_executable : str
         Path to the `simple3DFilter` executable [1].
     snapshot_kind : str
-        Kind of the simulation snapshot. Currently only `gadget4` is supported.
+        Kind of the simulation snapshot.
 
     Returns
     -------
@@ -147,12 +147,24 @@ def main(snapshot_path, output_path, resolution, scratch_space, SPH_executable,
     ----------
     [1] https://bitbucket.org/glavaux/cosmotool/src/master/sample/simple3DFilter.cpp  # noqa
     """
-    if snapshot_kind != "gadget4":
+    # First get the temporary file path.
-        raise NotImplementedError("Only GADGET HDF5 snapshots are supported.")
+    if snapshot_kind == "gadget4":
+        temporary_output_path = join(
+            scratch_space, generate_unique_id(snapshot_path))
+    elif snapshot_kind == "ramses":
+        temporary_output_path = snapshot_path
+    else:
+        raise NotImplementedError("Only GADGET HDF5 or preprocessed RAMSES "
+                                  "snapshots are supported.")
 
+    if not temporary_output_path.endswith(".hdf5"):
+        raise RuntimeError("Temporary output path must end with `.hdf5`.")
+
+    # Print the job information.
     print("---------- SPH Density & Velocity Field Job Information ----------")
     print(f"Snapshot path:     {snapshot_path}")
     print(f"Output path:       {output_path}")
+    print(f"Temporary path:    {temporary_output_path}")
     print(f"Resolution:        {resolution}")
     print(f"Scratch space:     {scratch_space}")
     print(f"SPH executable:    {SPH_executable}")
@@ -160,24 +172,28 @@ def main(snapshot_path, output_path, resolution, scratch_space, SPH_executable,
     print("------------------------------------------------------------------")
     print(flush=True)
 
-    temporary_output_path = join(
+    # Prepare or read-off the temporary snapshot file.
-        scratch_space, generate_unique_id(snapshot_path))
+    if snapshot_kind == "gadget4":
-
+        print(f"{now()}: preparing snapshot...", flush=True)
-    if not temporary_output_path.endswith(".hdf5"):
+        boxsize = prepare_gadget(snapshot_path, temporary_output_path)
-        raise RuntimeError("Temporary output path must end with `.hdf5`.")
+        print(f"{now()}: wrote temporary data to {temporary_output_path}.",
-
+              flush=True)
-    print(f"{now()}: preparing snapshot...", flush=True)
+    else:
-    boxsize = prepare_gadget(snapshot_path, temporary_output_path)
+        boxsize = 677.7  # Mpc/h
-    print(f"{now()}: wrote temporary data to {temporary_output_path}.",
+        print(f"{now()}: CAREFUL, forcefully setting the boxsize to {boxsize} Mpc / h.",  # noqa
-          flush=True)
+              flush=True)
 
+    # Run the SPH filter.
     run_sph_filter(temporary_output_path, output_path, boxsize, resolution,
                    SPH_executable)
-    print(f"{now()}: removing the temporary snapshot file.", flush=True)
+
-    try:
+    # Remove the temporary snapshot file if it was created.
-        remove(temporary_output_path)
+    if snapshot_kind == "gadget4":
-    except FileNotFoundError:
+        print(f"{now()}: removing the temporary snapshot file.", flush=True)
-        pass
+        try:
+            remove(temporary_output_path)
+        except FileNotFoundError:
+            pass
 
 
 if __name__ == "__main__":
@@ -193,7 +209,7 @@ if __name__ == "__main__":
     parser.add_argument("--SPH_executable", type=str, required=True,
                         help="Path to the `simple3DFilter` executable.")
     parser.add_argument("--snapshot_kind", type=str, required=True,
-                        choices=["gadget4"],
+                        choices=["gadget4", "ramses"],
                         help="Kind of the simulation snapshot.")
     args = parser.parse_args()
 

scripts_independent/run_field_sph_gadget.py

@@ -17,12 +17,27 @@ Script to write the SLURM submission script and submit it to the queue to
 calculate the SPH density & velocity field for GADGET.
 """
 from os import system
+from os.path import join
 
 
-def write_submit(chain_index, kind, resolution, nthreads):
+def write_submit(chain_index, kind, resolution, nthreads, snapshot_kind):
     if kind not in ["main", "random", "varysmall"]:
         raise RuntimeError(f"Unknown kind `{kind}`.")
 
+    basepath = "/cosma8/data/dp016/dc-stis1/"
+    if snapshot_kind == "gadget4":
+        snapshot_path = join(basepath, f"csiborg2_{kind}/chain_{chain_index}",
+                             "output/snapshot_099_full.hdf5")
+        output_path = join(
+            basepath,
+            f"csiborg2_{kind}/field/chain_{chain_index}_{resolution}.hdf5")
+    else:
+        chain_index = str(chain_index).zfill(5)
+        snapshot_path = join(basepath, "csiborg1_sph",
+                             f"ramses_{chain_index}.hdf5")
+        output_path = join(basepath, "csiborg1_sph",
+                           f"sph_ramses_{chain_index}_{resolution}.hdf5")
+
     txt = f"""#!/bin/sh
 
 #SBATCH --ntasks-per-node=1
@@ -53,12 +68,12 @@ source /cosma/home/dp016/dc-stis1/csiborgtools/venv_csiborgtools/bin/activate
 export OMP_NUM_THREADS={nthreads}
 export OMP_NESTED=true
 
-snapshot_path="/cosma8/data/dp016/dc-stis1/csiborg2_{kind}/chain_{chain_index}/output/snapshot_099_full.hdf5"
+snapshot_path={snapshot_path}
-output_path="/cosma8/data/dp016/dc-stis1/csiborg2_{kind}/field/chain_{chain_index}_{resolution}.hdf5"
+output_path={output_path}
 resolution={resolution}
 scratch_space="/snap8/scratch/dp016/dc-stis1/"
 SPH_executable="/cosma8/data/dp016/dc-stis1/cosmotool/bld2/sample/simple3DFilter"
-snapshot_kind="gadget4"
+snapshot_kind={snapshot_kind}
 
 python3 field_sph_gadget.py --snapshot_path $snapshot_path --output_path $output_path --resolution $resolution --scratch_space $scratch_space --SPH_executable $SPH_executable --snapshot_kind $snapshot_kind
 """
@@ -72,18 +87,25 @@ python3 field_sph_gadget.py --snapshot_path $snapshot_path --output_path $output
 
 
 if __name__ == "__main__":
+    snapshot_kind = "gadget4"
     # kind = "main"
     # chains = [15617, 15717, 15817, 15917, 16017, 16117, 16217, 16317, 16417, 16517, 16617, 16717, 16817, 16917, 17017, 17117, 17217, 17317, 17417]
 
     # kind = "varysmall"
     # chains = ["16417_001", "16417_025", "16417_050", "16417_075", "16417_100", "16417_125", "16417_150", "16417_175", "16417_200", "16417_225", "16417_250", "16417_275", "16417_300", "16417_325", "16417_350", "16417_375", "16417_400", "16417_425", "16417_450", "16417_475"]
 
-    kind = "random"
+    # kind = "random"
-    chains = [1, 25, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300, 325, 350, 375, 400, 425, 450, 475]
+    # chains = [1, 25, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300, 325, 350, 375, 400, 425, 450, 475]
+
+    # CSiBORG1 RAMSES
+    snapshot_kind = "ramses"
+    kind = "main"
+    chains = [7444]
 
     resolution = 1024
     nthreads = 32
 
     for chain_index in chains:
-        fname = write_submit(chain_index, kind, resolution, nthreads)
+        fname = write_submit(chain_index, kind, resolution, nthreads,
+                             snapshot_kind)
         system(f"sbatch {fname}")