commit 04cea3ab51deef1d328d4e7e23542f9d9f2b75f6 Author: glatterf42 Date: Wed Apr 13 13:15:11 2022 +0200 first commit diff --git a/__pycache__/directories.cpython-37.pyc b/__pycache__/directories.cpython-37.pyc new file mode 100644 index 0000000..fb6ade7 Binary files /dev/null and b/__pycache__/directories.cpython-37.pyc differ diff --git a/__pycache__/directories.cpython-38.pyc b/__pycache__/directories.cpython-38.pyc new file mode 100644 index 0000000..8b39588 Binary files /dev/null and b/__pycache__/directories.cpython-38.pyc differ diff --git a/__pycache__/write_spectra_jobs.cpython-37.pyc b/__pycache__/write_spectra_jobs.cpython-37.pyc new file mode 100644 index 0000000..e8a566c Binary files /dev/null and b/__pycache__/write_spectra_jobs.cpython-37.pyc differ diff --git a/__pycache__/write_spectra_jobs.cpython-38.pyc b/__pycache__/write_spectra_jobs.cpython-38.pyc new file mode 100644 index 0000000..9fca6ad Binary files /dev/null and b/__pycache__/write_spectra_jobs.cpython-38.pyc differ diff --git a/agora_plot_power_spectra.py b/agora_plot_power_spectra.py new file mode 100644 index 0000000..6959d09 --- /dev/null +++ b/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") + + + + diff --git a/directories.py b/directories.py new file mode 100644 index 0000000..cf7f880 --- /dev/null +++ b/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/') \ No newline at end of file diff --git a/plot_cross_spectra.py b/plot_cross_spectra.py new file mode 100644 index 0000000..b49be87 --- /dev/null +++ b/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") + + + + diff --git a/plot_power_spectra.py b/plot_power_spectra.py new file mode 100644 index 0000000..3a46db8 --- /dev/null +++ b/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") + + + + diff --git a/plot_spectra.py b/plot_spectra.py new file mode 100644 index 0000000..4c2166e --- /dev/null +++ b/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") + + + diff --git a/prepare_spectral_evaluation_agora.py b/prepare_spectral_evaluation_agora.py new file mode 100644 index 0000000..bd40582 --- /dev/null +++ b/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 + ) + diff --git a/prepare_spectral_evaluation_wavelets.py b/prepare_spectral_evaluation_wavelets.py new file mode 100644 index 0000000..2ad22b2 --- /dev/null +++ b/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 + ) + diff --git a/visualise_ics.py b/visualise_ics.py new file mode 100644 index 0000000..617ea5c --- /dev/null +++ b/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]) \ No newline at end of file diff --git a/write_spectra_jobs.py b/write_spectra_jobs.py new file mode 100644 index 0000000..4b3f6c2 --- /dev/null +++ b/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/" + ) \ No newline at end of file