From 0ee7c9591b6c2a5c238dbfe06af8efcb54fd7fcf Mon Sep 17 00:00:00 2001 From: thomashelfer Date: Mon, 28 Sep 2020 15:11:54 +0100 Subject: [PATCH 1/5] Added spherical integrator --- .../parallel_spherical_integrator.py | 160 ++++++++++++++++++ 1 file changed, 160 insertions(+) create mode 100644 YTAnalysisTools/parallel_spherical_integrator.py diff --git a/YTAnalysisTools/parallel_spherical_integrator.py b/YTAnalysisTools/parallel_spherical_integrator.py new file mode 100644 index 0000000..4dd8603 --- /dev/null +++ b/YTAnalysisTools/parallel_spherical_integrator.py @@ -0,0 +1,160 @@ +# parallel_Psi4.py +# James Widdicombe +# Last Updated 16/12/2019 +# Decomposition of Psi4 into spin weighted spherical harmonics +# l = 2,3,4 + +# Import the modules +import yt +import numpy as np +from numpy import pi +import time +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt +from matplotlib import rcParams + +start_time = time.time() + +# Matplotlib Settings +rcParams.update({"figure.autolayout": True}) +rcParams["axes.formatter.limits"] = [-3, 3] + +# Enable YT Parallelism +yt.enable_parallelism() + +# Script Parameters +extraction_radius = 60 # radius of extraction +RealPart = "chi" +# Imaginary part of the variable you want to integrate over +#ImagPart = "chi" +data_location = "../../ScalarField_*.3d.hdf5" # Data file location + +# Loading dataset +ts = yt.load(data_location) + +# Define Deltat for power output +time0 = ts[0].current_time +time1 = ts[1].current_time +DeltaT = float(time1 - time0) + +# Define an empty storage dictionary for collecting information +# in parallel through processing +storage = {} + +# Program Parameters +center = ts[0].domain_right_edge / 2.0 + +# Definitions for Quadrature scheme +N = 131 + +coord = np.loadtxt("PointDistFiles/lebedev/lebedev_%03d.txt" % N) +theta = coord[:, 1] * pi / 180 +phi = coord[:, 0] * pi / 180 + pi +w = coord[:, 2] + +phi_length = len(phi) + +# Iterate through dataset +for sto, i in ts.piter(storage=storage): + # Timings + L_start = time.time() + + # Initalising + Sr = 0 + 1j * 0 + # l = 2 + Weyl4_l2_m0 = 0 + 1j * 0 + + Weyl4data = [] + + # k is a counter + for k in range(phi_length): + + phi_var = phi[k] + theta_var = theta[k] + x1 = extraction_radius * np.cos(phi_var) * np.sin(theta_var) + float(center[0]) + y1 = extraction_radius * np.sin(phi_var) * np.sin(theta_var) + float(center[1]) + z1 = extraction_radius * np.cos(theta_var) + float(center[2]) + c = [x1, y1, z1] + ptn = i.point(c) + Real = float(ptn[RealPart][0]) + if ImagPart is not None: + Imag = float(ptn[ImagPart][0]) + else: + Imag = 0 + Integrand = Real + 1j * Imag + + + Weyl4_l2_m0 += ( + 4 * pi * w[k] * Integrand * extraction_radius ** 2 + ) + + # positive m + array = [ + i.current_time, + Weyl4_l2_m0, + time.time() - L_start, + ] + sto.result = array + sto.result_id = str(i) +if yt.is_root(): + + # Initialize Output arrays + # l = 2 + + timedata = [] + + # Diagnostics + loop_times = [] + + Integrated_values = [] + + # Swap from storage into arrays + for L in sorted(storage.items()): + timedata.append(L[1][0]) + Integrated_values.append(L[1][1]) + loop_times.append(L[1][2]) + + All_data = [] + All_data.append(Integrated_values) + + # Output Data + np.savetxt("time.out", timedata) + np.savetxt("loop_times.out", loop_times) + np.savetxt("speed_up.out", [sum(loop_times) / (time.time() - start_time)]) + np.savetxt("Integrated_values.out", Integrated_values) + + # Integrated Plotting + labels = [ + "integrated_evolution ", + ] + + # Calculate Retarded Time + time_retarded = [] + for i in timedata: + time_retarded.append(float(i) - extraction_radius) + np.savetxt("timeretarded.out", time_retarded) + + # Individual Modes + for i in range(0, len(labels)): + plt.figure(i) + plt.plot(time_retarded, np.real(All_data[i])) + plt.xlabel(r"$t_{ret}~[1/m]$") + plt.ylabel(r"$r\Psi_4$") + plt.grid() + plt.savefig(labels[i] + ".png", bbox_inches="tight") + plt.close() + # All Modes + plt.figure(len(labels), figsize=(10, 6)) + ax = plt.subplot(111) + for i in range(0, len(labels)): + ax.plot(time_retarded, np.real(All_data[i]), label=labels[i]) + box = ax.get_position() + ax.set_position([box.x0, box.y0, box.width * 0.8, box.height]) + plt.xlabel(r"$t_{ret}~[1/m_{a}]$") + plt.ylabel(r"$r\Psi_4$") + plt.grid() + plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0) + plt.savefig("All.png", bbox_inches="tight") + plt.close() From 0e6a9f392436f8ac0f9475f7557351b8b34f6bec Mon Sep 17 00:00:00 2001 From: thomashelfer Date: Mon, 28 Sep 2020 15:15:39 +0100 Subject: [PATCH 2/5] fixed minor bug --- YTAnalysisTools/parallel_spherical_integrator.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/YTAnalysisTools/parallel_spherical_integrator.py b/YTAnalysisTools/parallel_spherical_integrator.py index 4dd8603..a14dd9b 100644 --- a/YTAnalysisTools/parallel_spherical_integrator.py +++ b/YTAnalysisTools/parallel_spherical_integrator.py @@ -26,9 +26,10 @@ # Script Parameters extraction_radius = 60 # radius of extraction +# Name of variable RealPart = "chi" # Imaginary part of the variable you want to integrate over -#ImagPart = "chi" +ImagPart = None data_location = "../../ScalarField_*.3d.hdf5" # Data file location # Loading dataset From 4da7a028c47c6a7573769a162e8955895de77176 Mon Sep 17 00:00:00 2001 From: Thomas Helfer Date: Mon, 28 Sep 2020 16:23:02 +0200 Subject: [PATCH 3/5] autopep8 applied to spherial integrator --- .../parallel_spherical_integrator.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/YTAnalysisTools/parallel_spherical_integrator.py b/YTAnalysisTools/parallel_spherical_integrator.py index a14dd9b..65076fe 100644 --- a/YTAnalysisTools/parallel_spherical_integrator.py +++ b/YTAnalysisTools/parallel_spherical_integrator.py @@ -1,10 +1,11 @@ -# parallel_Psi4.py +# parallel_spherical_integrator.py # James Widdicombe -# Last Updated 16/12/2019 -# Decomposition of Psi4 into spin weighted spherical harmonics -# l = 2,3,4 +# Last Updated 28/09/2020 +# Integrate a field over a sphere # Import the modules +from matplotlib import rcParams +import matplotlib.pyplot as plt import yt import numpy as np from numpy import pi @@ -12,8 +13,6 @@ import matplotlib matplotlib.use("Agg") -import matplotlib.pyplot as plt -from matplotlib import rcParams start_time = time.time() @@ -74,8 +73,10 @@ phi_var = phi[k] theta_var = theta[k] - x1 = extraction_radius * np.cos(phi_var) * np.sin(theta_var) + float(center[0]) - y1 = extraction_radius * np.sin(phi_var) * np.sin(theta_var) + float(center[1]) + x1 = extraction_radius * \ + np.cos(phi_var) * np.sin(theta_var) + float(center[0]) + y1 = extraction_radius * \ + np.sin(phi_var) * np.sin(theta_var) + float(center[1]) z1 = extraction_radius * np.cos(theta_var) + float(center[2]) c = [x1, y1, z1] ptn = i.point(c) @@ -86,7 +87,6 @@ Imag = 0 Integrand = Real + 1j * Imag - Weyl4_l2_m0 += ( 4 * pi * w[k] * Integrand * extraction_radius ** 2 ) From 02fb4d6fb20b13599edc195c0d57c6800cfd077f Mon Sep 17 00:00:00 2001 From: thomashelfer Date: Fri, 2 Oct 2020 16:31:46 +0100 Subject: [PATCH 4/5] minor changes --- .../parallel_spherical_integrator.py | 17 ++--------------- 1 file changed, 2 insertions(+), 15 deletions(-) diff --git a/YTAnalysisTools/parallel_spherical_integrator.py b/YTAnalysisTools/parallel_spherical_integrator.py index 65076fe..18a332c 100644 --- a/YTAnalysisTools/parallel_spherical_integrator.py +++ b/YTAnalysisTools/parallel_spherical_integrator.py @@ -87,14 +87,14 @@ Imag = 0 Integrand = Real + 1j * Imag - Weyl4_l2_m0 += ( + Integral += ( 4 * pi * w[k] * Integrand * extraction_radius ** 2 ) # positive m array = [ i.current_time, - Weyl4_l2_m0, + Integral, time.time() - L_start, ] sto.result = array @@ -146,16 +146,3 @@ plt.grid() plt.savefig(labels[i] + ".png", bbox_inches="tight") plt.close() - # All Modes - plt.figure(len(labels), figsize=(10, 6)) - ax = plt.subplot(111) - for i in range(0, len(labels)): - ax.plot(time_retarded, np.real(All_data[i]), label=labels[i]) - box = ax.get_position() - ax.set_position([box.x0, box.y0, box.width * 0.8, box.height]) - plt.xlabel(r"$t_{ret}~[1/m_{a}]$") - plt.ylabel(r"$r\Psi_4$") - plt.grid() - plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0) - plt.savefig("All.png", bbox_inches="tight") - plt.close() From 8314ef731badee4b4a0bd25800734b8059687269 Mon Sep 17 00:00:00 2001 From: thomashelfer Date: Fri, 2 Oct 2020 16:37:37 +0100 Subject: [PATCH 5/5] renamed and fixed minor bug --- ...cal_integrator.py => parallel_integrator_over_sphere.py} | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) rename YTAnalysisTools/{parallel_spherical_integrator.py => parallel_integrator_over_sphere.py} (93%) diff --git a/YTAnalysisTools/parallel_spherical_integrator.py b/YTAnalysisTools/parallel_integrator_over_sphere.py similarity index 93% rename from YTAnalysisTools/parallel_spherical_integrator.py rename to YTAnalysisTools/parallel_integrator_over_sphere.py index 18a332c..66c03b2 100644 --- a/YTAnalysisTools/parallel_spherical_integrator.py +++ b/YTAnalysisTools/parallel_integrator_over_sphere.py @@ -61,10 +61,8 @@ # Timings L_start = time.time() - # Initalising - Sr = 0 + 1j * 0 - # l = 2 - Weyl4_l2_m0 = 0 + 1j * 0 + # Init + Integral = 0 + 1j * 0 Weyl4data = []