Separated collapse plot modelling and plotting functionality.

numerical/models/slot_h_collapse.py

 import os
 
-import matplotlib.pyplot as plt
-import math
 import numpy as np
 
 import numerical.bem as bem
 import numerical.util.gen_utils as gen
-import common.util.plotting_utils as pu
-import matplotlib.patches as patches
+from common.util.file_utils import lists_to_csv
+from common.util.plotting_utils import initialize_plt
 
-plt.rc('text', usetex=True)
-font = {'family': 'serif', 'size': 10, 'serif': ['cmr10']}
-plt.rc('font', **font)
-plt.rc('lines', linewidth=1, markersize=3)
-plt.rc('axes', linewidth=0.5)
-plt.rc('patch', linewidth=0.5)
+initialize_plt()
 
 Hs = np.linspace(1, 5, 5)
 W = 2
@@ -24,22 +17,8 @@ Y = 2
 normalize = True
 
 m_0 = 1
-n = 20000
+n = 2000
 
-fig_width = 5.31445
-fig = plt.figure(figsize=(fig_width, fig_width / 2))
-
-left_pos = 0.11
-right_pos = 0.98
-wspace = 0.3
-plt.subplots_adjust(top=0.95, bottom=0.3, left=left_pos, right=right_pos, wspace=wspace)
-fig.patch.set_facecolor('white')
-ax = fig.add_subplot(121)
-ax.locator_params(nbins=6)
-norm_ax = fig.add_subplot(122)
-norm_ax.locator_params(nbins=6)
-norm_ax.set_xlim(-5.5, 5.5)
-norm_ax.set_ylim(-1.1, 1.1)
 xs = Xs / (0.5 * W)
 
 for i, H in enumerate(Hs):
@@ -55,54 +34,9 @@ for i, H in enumerate(Hs):
     points[:, 1] = Y
     points[:, 2] = 0
     vels = bem.get_jet_dirs(points, centroids, normals, areas, m_0, R_inv, verbose=True)
-    theta_js = np.arctan2(vels[:, 1], vels[:, 0]) + 0.5 * np.pi
-
-    x_to_plot = xs
-    theta_j_to_plot = theta_js
-
-    theta_j_star, x_star = sorted(zip(theta_js, xs), key=lambda k: k[0])[-1]
-    norm_x_to_plot = np.divide(xs, x_star)
-    norm_theta_j_to_plot = np.divide(theta_js, theta_j_star)
-
-    label_frac = "h"
-    if i == 0:
-        ax.plot(x_to_plot, theta_j_to_plot, label=f"${label_frac} = {H / W:.2f}$")
-        norm_ax.plot(norm_x_to_plot, norm_theta_j_to_plot, label=f"${label_frac} = {H / W:.2f}$")
-    else:
-        ax.plot(x_to_plot, theta_j_to_plot, label=f"${label_frac} = {H / W:.2f}$", linestyle="--",
-                dashes=(i, 2 * i))
-        norm_ax.plot(norm_x_to_plot, norm_theta_j_to_plot, label=f"${label_frac} = {H / W:.2f}$", linestyle="--",
-                     dashes=(i, 2 * i))
-
-# TODO: Convert to data generator + data plotter
-
-label_pad = 0
-norm_ax.set_xlabel("$\\hat{x}$", labelpad=label_pad)
-norm_ax.set_ylabel("$\\hat{\\theta}$", labelpad=label_pad)
-ax.set_xlabel("$x$", labelpad=label_pad)
-ax.set_ylabel("$\\theta_j$", labelpad=label_pad)
-ax.axvline(x=-1, linestyle='--', color='gray')
-ax.axvline(x=1, linestyle='--', color='gray')
-# ax.legend()
-
-ax.legend(bbox_to_anchor=(0, -0.34, 2.3, .05), loc=10, ncol=len(Hs), mode="expand",
-          borderaxespad=0,
-          fancybox=False, edgecolor='k', shadow=False, handlelength=1.5, handletextpad=0.5)
-ax.annotate('($a$)', xy=(0, 0), xytext=(0.115, 0.89), textcoords='figure fraction', horizontalalignment='left',
-            verticalalignment='bottom')
-norm_ax.annotate('($b$)', xy=(0, 0), xytext=(0.61, 0.89), textcoords='figure fraction',
-                 horizontalalignment='left', verticalalignment='bottom')
-# ymin, ymax = ax.get_ylim()
-# ax.set_yticks(np.round(np.linspace(ymin, ymax, 5), 2))
-#
-# xmin, xmax = ax.get_xlim()
-# ax.set_xticks(np.round(np.linspace(xmin, xmax, 5), 2))
-#
-# ymin, ymax = norm_ax.get_ylim()
-# norm_ax.set_yticks(np.round(np.linspace(ymin, ymax, 5), 2))
-#
-# xmin, xmax = norm_ax.get_xlim()
-# norm_ax.set_xticks(np.round(np.linspace(xmin, xmax, 5), 2))
+    thetas = np.arctan2(vels[:, 1], vels[:, 0]) + 0.5 * np.pi
 
-plt.savefig('model_outputs/h_sweep_plot.pdf')
-plt.show()
+    f_dir = "model_outputs/slot_h_collapse/"
+    f_name = f"h_collapse_n{n}_H{H}.csv"
+    if not os.path.exists(f_dir + f_name):
+        lists_to_csv(f_dir, f_name, [xs, thetas], ["x", "theta"])

numerical/models/slot_q_collapse.py → numerical/models/slot_h_collapse_plot.py

@@ -6,25 +6,21 @@ import numpy as np
 
 import numerical.bem as bem
 import numerical.util.gen_utils as gen
-import common.util.plotting_utils as pu
+from common.util.plotting_utils import initialize_plt
+from common.util.file_utils import csv_to_lists
 import matplotlib.patches as patches
 
-plt.rc('text', usetex=True)
-font = {'family': 'serif', 'size': 10, 'serif': ['cmr10']}
-plt.rc('font', **font)
-plt.rc('lines', linewidth=1, markersize=3)
-plt.rc('axes', linewidth=0.5)
-plt.rc('patch', linewidth=0.5)
+initialize_plt()
 
-H = 2
+Hs = np.linspace(1, 5, 5)
 W = 2
 Xs = np.linspace(-3 * W, 3 * W, 300)
-Ys = np.linspace(1, 5, 5)
+Y = 2
 
 normalize = True
 
 m_0 = 1
-n = 20000
+n = 2000
 
 fig_width = 5.31445
 fig = plt.figure(figsize=(fig_width, fig_width / 2))
@@ -42,49 +38,36 @@ norm_ax.set_xlim(-5.5, 5.5)
 norm_ax.set_ylim(-1.1, 1.1)
 xs = Xs / (0.5 * W)
 
-centroids, normals, areas = gen.gen_varied_slot(n=n, H=H, W=W, length=50, depth=50, w_thresh=12, density_ratio=0.25)
-print("Requested n = {0}, using n = {1}.".format(n, len(centroids)))
-print(np.mean(centroids, 0))
-R_matrix = bem.get_R_matrix(centroids, normals, areas, dtype=np.float32)
-R_inv = np.linalg.inv(R_matrix)
-for i, Y in enumerate(Ys):
-    print(f"Testing Y = {Y}")
-    points = np.empty((len(Xs), 3))
-    points[:, 0] = Xs
-    points[:, 1] = Y
-    points[:, 2] = 0
-    vels = bem.get_jet_dirs(points, centroids, normals, areas, m_0, R_inv, verbose=True)
-    theta_js = np.arctan2(vels[:, 1], vels[:, 0]) + 0.5 * np.pi
+for i, H in enumerate(Hs):
+    print(f"Reading H = {H}")
+    f_dir = "model_outputs/slot_h_collapse/"
+    f_name = f"h_collapse_n{n}_H{H}.csv"
+    xs, thetas = csv_to_lists(f_dir, f_name, True)
 
-    x_to_plot = xs
-    theta_j_to_plot = theta_js
-
-    theta_j_star, x_star = sorted(zip(theta_js, xs), key=lambda k: k[0])[-1]
+    theta_star, x_star = sorted(zip(thetas, xs), key=lambda k: k[0])[-1]
     norm_x_to_plot = np.divide(xs, x_star)
-    norm_theta_j_to_plot = np.divide(theta_js, theta_j_star)
+    norm_theta_to_plot = np.divide(thetas, theta_star)
 
-    label_frac = "y"
+    label_frac = "h"
     if i == 0:
-        ax.plot(x_to_plot, theta_j_to_plot, label=f"${label_frac} = {Y / W:.2f}$")
-        norm_ax.plot(norm_x_to_plot, norm_theta_j_to_plot, label=f"${label_frac} = {Y / W:.2f}$")
+        ax.plot(xs, thetas, label=f"${label_frac} = {H / W:.2f}$")
+        norm_ax.plot(norm_x_to_plot, norm_theta_to_plot, label=f"${label_frac} = {H / W:.2f}$")
     else:
-        ax.plot(x_to_plot, theta_j_to_plot, label=f"${label_frac} = {Y / W:.2f}$", linestyle="--",
+        ax.plot(xs, thetas, label=f"${label_frac} = {H / W:.2f}$", linestyle="--",
                 dashes=(i, 2 * i))
-        norm_ax.plot(norm_x_to_plot, norm_theta_j_to_plot, label=f"${label_frac} = {Y / W:.2f}$", linestyle="--",
+        norm_ax.plot(norm_x_to_plot, norm_theta_to_plot, label=f"${label_frac} = {H / W:.2f}$", linestyle="--",
                      dashes=(i, 2 * i))
 
-# TODO: Convert to data generator + data plotter
-
 label_pad = 0
 norm_ax.set_xlabel("$\\hat{x}$", labelpad=label_pad)
 norm_ax.set_ylabel("$\\hat{\\theta}$", labelpad=label_pad)
 ax.set_xlabel("$x$", labelpad=label_pad)
-ax.set_ylabel("$\\theta_j$", labelpad=label_pad)
+ax.set_ylabel("$\\theta$ (rad)", labelpad=label_pad)
 ax.axvline(x=-1, linestyle='--', color='gray')
 ax.axvline(x=1, linestyle='--', color='gray')
 # ax.legend()
 
-ax.legend(bbox_to_anchor=(0, -0.34, 2.3, .05), loc=10, ncol=len(Ys), mode="expand",
+ax.legend(bbox_to_anchor=(0, -0.34, 2.3, .05), loc=10, ncol=len(Hs), mode="expand",
           borderaxespad=0,
           fancybox=False, edgecolor='k', shadow=False, handlelength=1.5, handletextpad=0.5)
 ax.annotate('($a$)', xy=(0, 0), xytext=(0.115, 0.89), textcoords='figure fraction', horizontalalignment='left',
@@ -103,5 +86,5 @@ norm_ax.annotate('($b$)', xy=(0, 0), xytext=(0.61, 0.89), textcoords='figure fra
 # xmin, xmax = norm_ax.get_xlim()
 # norm_ax.set_xticks(np.round(np.linspace(xmin, xmax, 5), 2))
 
-plt.savefig('model_outputs/q_sweep_plot.pdf')
+plt.savefig('model_outputs/h_sweep_plot.pdf')
 plt.show()

numerical/models/slot_y_collapse.py

+import os
+
+import numpy as np
+
+import numerical.bem as bem
+import numerical.util.gen_utils as gen
+from common.util.file_utils import lists_to_csv
+
+H = 2
+W = 2
+Xs = np.linspace(-3 * W, 3 * W, 300)
+Ys = np.linspace(1, 5, 5)
+
+m_0 = 1
+n = 20000
+
+xs = Xs / (0.5 * W)
+
+centroids, normals, areas = gen.gen_varied_slot(n=n, H=H, W=W, length=50, depth=50, w_thresh=12, density_ratio=0.25)
+print("Requested n = {0}, using n = {1}.".format(n, len(centroids)))
+print(np.mean(centroids, 0))
+R_matrix = bem.get_R_matrix(centroids, normals, areas, dtype=np.float32)
+R_inv = np.linalg.inv(R_matrix)
+for i, Y in enumerate(Ys):
+    print(f"Testing Y = {Y}")
+    points = np.empty((len(Xs), 3))
+    points[:, 0] = Xs
+    points[:, 1] = Y
+    points[:, 2] = 0
+    vels = bem.get_jet_dirs(points, centroids, normals, areas, m_0, R_inv, verbose=True)
+    thetas = np.arctan2(vels[:, 1], vels[:, 0]) + 0.5 * np.pi
+
+    x_to_plot = xs
+    theta_to_plot = thetas
+
+    theta_star, x_star = sorted(zip(thetas, xs), key=lambda k: k[0])[-1]
+    norm_x_to_plot = np.divide(xs, x_star)
+    norm_theta_to_plot = np.divide(thetas, theta_star)
+
+    f_dir = "model_outputs/slot_y_collapse/"
+    f_name = f"y_collapse_n{n}_Y{Y}.csv"
+    if not os.path.exists(f_dir + f_name):
+        lists_to_csv(f_dir, f_name, [xs, thetas], ["x", "theta"])

numerical/models/slot_y_collapse_plot.py

+import os
+
+import matplotlib.pyplot as plt
+import math
+import numpy as np
+
+import numerical.bem as bem
+import numerical.util.gen_utils as gen
+from common.util.plotting_utils import initialize_plt
+from common.util.file_utils import csv_to_lists
+import matplotlib.patches as patches
+
+initialize_plt()
+
+H = 2
+W = 2
+Xs = np.linspace(-3 * W, 3 * W, 300)
+Ys = np.linspace(1, 5, 5)
+
+m_0 = 1
+n = 2000
+
+fig_width = 5.31445
+fig = plt.figure(figsize=(fig_width, fig_width / 2))
+
+left_pos = 0.11
+right_pos = 0.98
+wspace = 0.3
+plt.subplots_adjust(top=0.95, bottom=0.3, left=left_pos, right=right_pos, wspace=wspace)
+fig.patch.set_facecolor('white')
+ax = fig.add_subplot(121)
+ax.locator_params(nbins=6)
+norm_ax = fig.add_subplot(122)
+norm_ax.locator_params(nbins=6)
+norm_ax.set_xlim(-5.5, 5.5)
+norm_ax.set_ylim(-1.1, 1.1)
+xs = Xs / (0.5 * W)
+
+for i, Y in enumerate(Ys):
+    print(f"Reading Y = {Y}")
+    f_dir = "model_outputs/slot_y_collapse/"
+    f_name = f"y_collapse_n{n}_Y{Y}.csv"
+    xs, thetas = csv_to_lists(f_dir, f_name, True)
+
+    theta_star, x_star = sorted(zip(thetas, xs), key=lambda k: k[0])[-1]
+    norm_x_to_plot = np.divide(xs, x_star)
+    norm_theta_to_plot = np.divide(thetas, theta_star)
+
+    label_frac = "y"
+    if i == 0:
+        ax.plot(xs, thetas, label=f"${label_frac} = {Y / W:.2f}$")
+        norm_ax.plot(norm_x_to_plot, norm_theta_to_plot, label=f"${label_frac} = {Y / W:.2f}$")
+    else:
+        ax.plot(xs, thetas, label=f"${label_frac} = {Y / W:.2f}$", linestyle="--",
+                dashes=(i, 2 * i))
+        norm_ax.plot(norm_x_to_plot, norm_theta_to_plot, label=f"${label_frac} = {Y / W:.2f}$", linestyle="--",
+                     dashes=(i, 2 * i))
+
+label_pad = 0
+norm_ax.set_xlabel("$\\hat{x}$", labelpad=label_pad)
+norm_ax.set_ylabel("$\\hat{\\theta}$", labelpad=label_pad)
+ax.set_xlabel("$x$", labelpad=label_pad)
+ax.set_ylabel("$\\theta$ (rad)", labelpad=label_pad)
+ax.axvline(x=-1, linestyle='--', color='gray')
+ax.axvline(x=1, linestyle='--', color='gray')
+# ax.legend()
+
+ax.legend(bbox_to_anchor=(0, -0.34, 2.3, .05), loc=10, ncol=len(Ys), mode="expand",
+          borderaxespad=0,
+          fancybox=False, edgecolor='k', shadow=False, handlelength=1.5, handletextpad=0.5)
+ax.annotate('($a$)', xy=(0, 0), xytext=(0.05, 0.95), textcoords='axes fraction', horizontalalignment='left',
+            verticalalignment='top')
+norm_ax.annotate('($b$)', xy=(0, 0), xytext=(0.05, 0.95), textcoords='axes fraction',
+                 horizontalalignment='left', verticalalignment='top')
+# ymin, ymax = ax.get_ylim()
+# ax.set_yticks(np.round(np.linspace(ymin, ymax, 5), 2))
+#
+# xmin, xmax = ax.get_xlim()
+# ax.set_xticks(np.round(np.linspace(xmin, xmax, 5), 2))
+#
+# ymin, ymax = norm_ax.get_ylim()
+# norm_ax.set_yticks(np.round(np.linspace(ymin, ymax, 5), 2))
+#
+# xmin, xmax = norm_ax.get_xlim()
+# norm_ax.set_xticks(np.round(np.linspace(xmin, xmax, 5), 2))
+
+plt.show()