first commit

agora_plot_power_spectra.py

@@ -0,0 +1,88 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Mar 16 15:12:21 2022
+
+@author: ben
+"""
+
+import pandas as pd
+import matplotlib.pyplot as plt
+from pathlib import Path
+import scipy.special as sf
+import numpy as np
+
+def Dplus( lambda0, a ):
+	return a * np.sqrt(1.0+lambda0*a**3) * sf.hyp2f1(3.0/2.0, 5.0/6.0, 11.0/6.0, -lambda0 * a**3 )
+
+basedir = Path("/home/ben/sims/swiftsim/examples/agora/spectra/")
+
+#choose waveform and Lbox:
+Lbox = 100.0  #only option as of now
+Nres1 = 128
+# Nres2 = 256
+k0 = 2 * 3.14159265358979323846264338327950 / Lbox
+knyquist1 = Nres1 * k0
+# knyquist2 = Nres2 * k0
+times = [0.166666, 0.333333, 0.5, 0.666666, 1.0]
+# scale_factor = 0       # give index of a list above
+
+Omega_m  = 0.272
+Omega_L  = 0.728
+lambda_0 = Omega_L / Omega_m
+
+#find columns in file manually
+#is k really in Mpc? Swift doesn't use /h internally at least.
+columns = ["k [Mpc]", "Pcross", "P1", "err. P1", "P2", "err. P2", "P2-1", "err. P2-1", "modes in bin"]
+
+zstart = 100
+a_ics = 1 / (1 + zstart)
+filename_ics = basedir / 'agora_ics_cross_spectrum'
+df_ics = pd.read_csv(f'{filename_ics}.txt', sep=' ', skipinitialspace=True, header=None, names=columns, skiprows=1)
+#only consider rows above resolution limit
+df_ics = df_ics[df_ics['k [Mpc]'] >= k0]
+p1_ics = df_ics['P1']
+D_squared_ics = Dplus(lambda0=lambda_0, a=a_ics) ** 2
+p1_ics_noramlised = p1_ics / D_squared_ics
+
+for scale_factor in range(len(times)):
+    filename = basedir / f"agora_a{scale_factor}_cross_spectrum"
+# filename = basedir / f"{waveform}_{Lbox:.0f}/{waveform}_{Lbox:.0f}_ics_vsc_cross_spectrum" # for ICs
+# savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/power_{waveform}_{Lbox:.0f}_ics_vsc") # for ICs
+# plt.title(f"Power Spectra {waveform} L={Lbox:.0f} a=0.02 vsc") # for ICs
+
+
+    df = pd.read_csv(f"{filename}.txt", sep=" ", skipinitialspace=True, header=None, names=columns, skiprows=1)
+
+    #only consider rows above resolution limit
+    df = df[df["k [Mpc]"] >= k0]
+
+    k = df["k [Mpc]"]
+    p1 = df["P1"]
+    p1_error = df["err. P1"]
+    # p2 = df["P2"]
+    # p2_error = df["err. P2"]
+    # pcross = df["Pcross"]
+
+    D_squared = Dplus(lambda0=lambda_0, a=times[scale_factor]) ** 2
+    p1_normalised = p1 / D_squared
+    p1_ic_normalised = p1_normalised / p1_ics_noramlised
+
+    # Plot the power spectra:
+    plt.loglog(k, p1_ic_normalised, label=f"{times[scale_factor]}")
+    # plt.loglog(k, p2, label="P2")
+
+plt.title(f"Power Spectra Agora 128")
+savedir = Path(f"/home/ben/Pictures/swift/agora/spectra/power_spectra")
+
+plt.xlabel("k [Mpc]")
+plt.ylabel("P")
+plt.vlines(knyquist1, ymin=min(p1), ymax=max(p1), color="black", linestyles="dashed", label=f"k_ny {Nres1}")
+# plt.vlines(knyquist2, ymin=min(p2), ymax=max(p2), color="black", linestyles="dashed", label=f"{Nres2}")
+plt.legend()
+
+plt.savefig(f"{savedir}.png")
+
+
+
+

directories.py

@@ -0,0 +1,15 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 15 12:11:41 2022
+
+@author: ben
+"""
+
+from pathlib import Path
+
+spectra_basedir = Path("/home/ben/sims/spectra/")
+
+monofonic_tests_basedir = Path("/home/ben/sims/swift/monofonic_tests/")
+
+agora_test_basedir = Path('/home/ben/sims/swiftsim/examples/agora/')

plot_cross_spectra.py

@@ -0,0 +1,67 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Mar 16 15:12:21 2022
+
+@author: ben
+"""
+
+import pandas as pd
+import matplotlib.pyplot as plt
+from pathlib import Path
+
+basedir = Path("/home/ben/sims/swift/monofonic_tests/spectra/")
+
+#choose Nres and Lbox:
+waveforms = ['DB2', "DB4", "DB8", "shannon"]  #DB2, DB4, DB8, shannon are all we have right now
+Lbox = 100.0  #only option as of now
+Nres = 256      #128 and 256 exist for now
+k0 = 2 * 3.14159265358979323846264338327950 / Lbox
+knyquist = Nres * k0
+a = [0.166666, 0.333333, 0.5, 0.666666, 1.0]
+scale_factor = 4       # give index of a list above
+
+
+for wave in waveforms:
+    filename = basedir / f"{wave}_{Lbox:.0f}/{wave}_{Lbox:.0f}_a{scale_factor}_cross_spectrum"
+    # filename = basedir / f"{wave}_{Lbox:.0f}/{wave}_{Lbox:.0f}_ics_vsc_cross_spectrum" # for ICs
+    
+    #find columns in file manually
+    #is k really in Mpc? Swift doesn't use /h internally at least.
+    columns = ["k [Mpc]", "Pcross", "P1", "err. P1", "P2", "err. P2", "P2-1", "err. P2-1", "modes in bin"]
+    
+    df = pd.read_csv(f"{filename}.txt", sep=" ", skipinitialspace=True, header=None, names=columns, skiprows=1)
+    
+    #only consider rows above resolution limit
+    df = df[df["k [Mpc]"] >= k0]
+    
+    k = df["k [Mpc]"]
+    p1 = df["P1"]
+    p1_error = df["err. P1"]
+    p2 = df["P2"]
+    p2_error = df["err. P2"]
+    pcross = df["Pcross"]
+    
+    
+    # Plot the Cross Correlation:
+    plt.plot(k, pcross, label=f"{wave}")
+
+# savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/cross_{Nres}_{Lbox:.0f}_ics_vsc") # for ICs
+# plt.title(f"Cross correlation N={Nres} L={Lbox:.0f} a=0.02 vsc") # for ICs
+
+savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/cross_{Nres}_{Lbox:.0f}_a{scale_factor}")
+plt.title(f"Cross correlation N={Nres} L={Lbox:.0f} a={a[scale_factor]}")
+
+plt.xscale("log")
+plt.xlabel("k [Mpc]")
+plt.ylabel("C = Pcross")
+plt.ylim(0.8, 1.0)
+plt.xlim(k[0], knyquist)
+# plt.vlines(knyquist, ymin=min(p1), ymax=max(p1), color="black", linestyles="dashed", label=f"{Nres}")
+plt.legend()
+
+plt.savefig(f"{savedir}.png")
+
+
+
+

plot_power_spectra.py

@@ -0,0 +1,65 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Mar 16 15:12:21 2022
+
+@author: ben
+"""
+
+import pandas as pd
+import matplotlib.pyplot as plt
+from pathlib import Path
+
+basedir = Path("/home/ben/sims/swift/monofonic_tests/spectra/")
+
+#choose waveform and Lbox:
+waveform = "shannon"  #DB2, DB4, DB8 or shannon
+Lbox = 100.0  #only option as of now
+Nres1 = 128
+Nres2 = 256
+k0 = 2 * 3.14159265358979323846264338327950 / Lbox
+knyquist1 = Nres1 * k0
+knyquist2 = Nres2 * k0
+a = [0.166666, 0.333333, 0.5, 0.666666, 1.0]
+scale_factor = 4       # give index of a list above
+
+filename = basedir / f"{waveform}_{Lbox:.0f}/{waveform}_{Lbox:.0f}_a{scale_factor}_cross_spectrum"
+# filename = basedir / f"{waveform}_{Lbox:.0f}/{waveform}_{Lbox:.0f}_ics_vsc_cross_spectrum" # for ICs
+# savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/power_{waveform}_{Lbox:.0f}_ics_vsc") # for ICs
+# plt.title(f"Power Spectra {waveform} L={Lbox:.0f} a=0.02 vsc") # for ICs
+
+#find columns in file manually
+#is k really in Mpc? Swift doesn't use /h internally at least.
+columns = ["k [Mpc]", "Pcross", "P1", "err. P1", "P2", "err. P2", "P2-1", "err. P2-1", "modes in bin"]
+
+df = pd.read_csv(f"{filename}.txt", sep=" ", skipinitialspace=True, header=None, names=columns, skiprows=1)
+
+#only consider rows above resolution limit
+df = df[df["k [Mpc]"] >= k0]
+
+k = df["k [Mpc]"]
+p1 = df["P1"]
+p1_error = df["err. P1"]
+p2 = df["P2"]
+p2_error = df["err. P2"]
+pcross = df["Pcross"]
+
+
+# Plot the power spectra:
+plt.loglog(k, p1, label="P1")
+plt.loglog(k, p2, label="P2")
+
+plt.title(f"Power Spectra {waveform} L={Lbox:.0f} a={a[scale_factor]}")
+savedir = Path(f"/home/ben/Pictures/swift/monofonic_tests/spectra/power_{waveform}_{Lbox:.0f}_a{scale_factor}")
+
+plt.xlabel("k [Mpc]")
+plt.ylabel("P")
+plt.vlines(knyquist1, ymin=min(p1), ymax=max(p1), color="black", linestyles="dashed", label=f"{Nres1}")
+plt.vlines(knyquist2, ymin=min(p2), ymax=max(p2), color="black", linestyles="dashed", label=f"{Nres2}")
+plt.legend()
+
+plt.savefig(f"{savedir}.png")
+
+
+
+

plot_spectra.py

@@ -0,0 +1,60 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Mar 16 15:12:21 2022
+
+@author: ben
+"""
+
+import pandas as pd
+import matplotlib.pyplot as plt
+from pathlib import Path
+
+basedir = Path("/home/ben/sims/swift/monofonic_tests/spectra/")
+
+#choose waveform and Lbox:
+waveform = "DB8"  #DB4, DB8 or shannon
+Lbox = 100.0  #only option as of now
+Nres1 = 128
+Nres2 = 256
+k0 = 2 * 3.14159265358979323846264338327950 / Lbox
+knyquist1 = Nres1 * k0
+knyquist2 = Nres2 * k0
+
+filename = basedir / f"{waveform}_{Lbox:.0f}/{waveform}_{Lbox:.0f}_cross_spectrum"
+
+#find columns in file manually
+#is k really in Mpc? Swift doesn't use /h internally at least.
+columns = ["k [Mpc]", "Pcross", "P1", "err. P1", "P2", "err. P2", "P2-1", "err. P2-1", "modes in bin"]
+
+df = pd.read_csv(f"{filename}.txt", sep=" ", skipinitialspace=True, header=None, names=columns, skiprows=1)
+
+#only consider rows above resolution limit
+df = df[df["k [Mpc]"] >= k0]
+
+k = df["k [Mpc]"]
+p1 = df["P1"]
+p1_error = df["err. P1"]
+p2 = df["P2"]
+p2_error = df["err. P2"]
+pcross = df["Pcross"]
+
+
+# Plot the power spectra:
+plt.loglog(k, p1)
+plt.loglog(k, p2)
+plt.title(f"Power Spectra {waveform} {Lbox:.0f}")
+plt.xlabel("k [Mpc]")
+plt.ylabel("P")
+plt.vlines(knyquist1, ymin=min(p1), ymax=max(p1), color="black", linestyles="dashed", label=f"{Nres1}")
+plt.vlines(knyquist2, ymin=min(p2), ymax=max(p2), color="black", linestyles="dashed", label=f"{Nres2}")
+plt.legend()
+
+# Plot the cross correlation:
+# plt.plot(k, pcross)
+# plt.title(f"Cross correlation {waveform} {Lbox:.0f}")
+# plt.xlabel("k [Mpc]")
+# plt.ylabel("C = Pcross")
+
+
+

prepare_spectral_evaluation_agora.py

@@ -0,0 +1,34 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 15 12:00:49 2022
+
+@author: ben
+"""
+
+from sys import argv
+from pathlib import Path
+from write_spectra_jobs import write_spectra_jobs_agora
+from directories import monofonic_tests_basedir, agora_test_basedir
+
+def main(scale_factor: int, Nres1: int, Nres2: int):
+    a = [0.166666, 0.333333, 0.5, 0.666666, 1.0]
+    print(f"Chose scale factor a = {a[scale_factor]} (index {scale_factor} / {len(a) - 1})")
+    
+    write_spectra_jobs_agora(scale_factor, Nres1, Nres2, 
+        output_basedir=agora_test_basedir / "spectra/")
+    
+    
+    
+    
+if __name__ == "__main__":    
+    scale_factor = int(argv[1])
+    Nres1 = int(argv[2])
+    Nres2 = int(argv[3])
+    
+    main(
+        scale_factor=scale_factor,
+        Nres1=Nres1,
+        Nres2=Nres2
+        )
+

prepare_spectral_evaluation_wavelets.py

@@ -0,0 +1,45 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 15 12:00:49 2022
+
+@author: ben
+"""
+
+from sys import argv
+from pathlib import Path
+from write_spectra_jobs import write_spectra_jobs
+from directories import monofonic_tests_basedir
+
+def main(scale_factor: int, Nres1: int, Nres2: int, Lbox: float, waveforms: list):
+    a = [0.166666, 0.333333, 0.5, 0.666666, 1.0]
+    print(f"Chose scale factor a = {a[scale_factor]} (index {scale_factor} / {len(a) - 1})")
+    
+    for wave in waveforms:
+        write_spectra_jobs(scale_factor, Nres1, Nres2, Lbox, wave, 
+            output_basedir=monofonic_tests_basedir / "spectra/")
+    
+    
+    
+    
+if __name__ == "__main__":
+    if argv[5] == "all":
+        waveforms = ["DB2", "DB4", "DB6", "DB8", "DB10", "shannon"]
+    else:
+        waveforms = argv[5:len(argv)]
+        
+    assert len(waveforms) <= 6
+    
+    scale_factor = int(argv[1])
+    Nres1 = int(argv[2])
+    Nres2 = int(argv[3])
+    Lbox = float(argv[4])
+    
+    main(
+        scale_factor=scale_factor,
+        Nres1=Nres1,
+        Nres2=Nres2,
+        Lbox=Lbox,
+        waveforms=waveforms
+        )
+

visualise_ics.py

@@ -0,0 +1,26 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 21 10:13:36 2022
+
+@author: ben
+"""
+
+import h5py
+from pathlib import Path
+
+directory = Path(r"/home/ben/sims/swift/monofonic_tests/DB4_256_100/")
+file = h5py.File(directory / "ics_DB4_256_100_nan_test_40tasks.hdf5", "r") 
+
+# directory = Path(r"/home/ben/monofonic-experimental/")
+# file = h5py.File(directory / "ics_DB4_256_100.hdf5", "r")
+
+for key in file.keys():
+    print(key) #for finding all header entries, which are:
+    
+Header = file["Header"]
+ICs_parameters = file["ICs_parameters"]
+PartType1 = file["PartType1"]
+Units = file["Units"]
+
+print(PartType1['Coordinates'][0:10])

write_spectra_jobs.py

@@ -0,0 +1,87 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Mar 15 12:04:15 2022
+
+@author: ben
+
+This currently only evaluates the last snapshot.
+"""
+
+from sys import argv
+from pathlib import Path
+from directories import spectra_basedir, monofonic_tests_basedir, agora_test_basedir
+import os 
+import stat
+
+def write_spectra_jobs(scale_factor: int, Nres1: int, Nres2: int, Lbox: float, form: str, output_basedir: Path):
+    Nres_max = max(Nres1, Nres2)
+    input_dir_1 = str(monofonic_tests_basedir) + f"/{form}_{Nres1}_{Lbox:.0f}"
+    input_dir_2 = str(monofonic_tests_basedir) + f"/{form}_{Nres2}_{Lbox:.0f}" 
+    output_dir = output_basedir / f"{form}_{Lbox:.0f}"
+    output_dir_string = str(output_dir)
+    if output_dir.exists():
+        print(output_dir, "already exists. Skipping...")
+    else:
+        print("creating", output_dir)
+        output_dir.mkdir()
+    
+    
+    script = f"""#!/bin/bash
+set -u
+
+{spectra_basedir}/build/spectra --format=3 --output={output_dir_string}/{form}_{Lbox:.0f}_a{scale_factor} --ngrid={2 * Nres_max} --input={input_dir_1}/output_000{scale_factor}.hdf5 --input={input_dir_2}/output_000{scale_factor}.hdf5
+    """
+    
+    filename = output_dir / "generate_final_spectra.sh"
+    with (filename).open("w") as f:
+        f.write(script)
+
+    permissions = os.stat(filename)
+    os.chmod(filename, permissions.st_mode | stat.S_IEXEC)
+
+
+
+def write_spectra_jobs_agora(scale_factor: int, Nres1: int, Nres2: int, output_basedir: Path):
+    Nres_max = max(Nres1, Nres2)
+    input_dir = str(agora_test_basedir) 
+    output_dir = output_basedir
+    output_dir_string = str(output_dir)
+    if output_dir.exists():
+        print(output_dir, "already exists. Skipping...")
+    else:
+        print("creating", output_dir)
+        output_dir.mkdir()
+    
+    script = f"""#!/bin/bash
+set -u
+
+{spectra_basedir}/build/spectra --format=3 --output={output_dir_string}/agora_a{scale_factor} --ngrid={2 * Nres_max} --input={input_dir}/output_000{scale_factor}.hdf5 --input={input_dir}/output_000{scale_factor}.hdf5
+    """
+    
+    filename = output_dir / "generate_final_spectra.sh"
+    with (filename).open("w") as f:
+        f.write(script)
+
+    permissions = os.stat(filename)
+    os.chmod(filename, permissions.st_mode | stat.S_IEXEC)
+    
+
+    
+if __name__ == "__main__": 
+    if argv[5] == "all":
+        waveforms = ["DB2", "DB4", "DB6", "DB8", "DB10", "shannon"]
+    else:
+        waveforms = argv[5:len(argv)]
+        
+    assert len(waveforms) <= 6
+    
+    for form in waveforms:
+        write_spectra_jobs(
+            scale_factor = int(argv[1]),
+            Nres1 = int(argv[2]),
+            Nres2 = int(argv[3]),
+            Lbox = float(argv[4]),
+            form = form,
+            output_dir = monofonic_tests_basedir / "spectra/"
+            )