Commit 8b9236da authored by Elijah Andrews's avatar Elijah Andrews

Reorganised numerical models.

parent 2d90512f
......@@ -12,7 +12,7 @@ from scipy import stats
from experimental.plotting.analyse_slot import SweepData, select_data_series
import experimental.util.analysis_utils as au
from numerical.models.slot_opt import find_slot_peak
from numerical.models.slot.slot_opt import find_slot_peak
import numerical.util.gen_utils as gen
import numerical.bem as bem
from common.util.plotting_utils import initialize_plt
......
......@@ -2,7 +2,7 @@ import numpy as np
import itertools
import matplotlib.pyplot as plt
from numerical.models.peters_corner import get_peters_corner_jet_dir
from numerical.models.corner.peters_corner import get_peters_corner_jet_dir
n = 100
xs = np.linspace(0.05, 1, n)
......
......@@ -31,7 +31,7 @@ fig.patch.set_facecolor('white')
ax = plt.gca()
ax.plot(ps, theta_js, label="Numeric")
ax.set_xlabel("$\\theta_b$", fontsize=18)
ax.set_ylabel("$\\theta_j$", fontsize=18)
ax.set_ylabel("$\\theta_j$ (rad)", fontsize=18)
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels, loc=4)
ax.axvline(0, linestyle='--', color='gray')
......
......@@ -23,7 +23,7 @@ def plot_peak_sweep(h_mat, y_mat, theta_star_mat, x_star_mat, theta_grad_mat,
theta_star_surf = ax.plot_surface(h_new_mat, y_new_mat, tsm_new, cmap=plt.cm.get_cmap('coolwarm'))
ax.set_xlabel("$h$")
ax.set_ylabel("$y$")
ax.set_zlabel("$\\theta_j^\\star$")
ax.set_zlabel("$\\theta^\\star$ (rad)")
theta_star_surf.set_clim(np.nanmin(theta_star_mat), np.nanmax(theta_star_mat))
theta_star_fig.colorbar(theta_star_surf)
......@@ -33,7 +33,7 @@ def plot_peak_sweep(h_mat, y_mat, theta_star_mat, x_star_mat, theta_grad_mat,
cmap=plt.cm.get_cmap('coolwarm'))
ax.set_xlabel("$log_{10}(h)$")
ax.set_ylabel("$log_{10}(y)$")
ax.set_zlabel("$log_{10}(\\theta_j^\\star)$")
ax.set_zlabel("$log_{10}(\\theta^\\star)$")
log_theta_star_surf.set_clim(np.nanmin(theta_star_mat), np.nanmax(theta_star_mat))
log_theta_star_fig.colorbar(log_theta_star_surf)
......@@ -56,7 +56,7 @@ def plot_peak_sweep(h_mat, y_mat, theta_star_mat, x_star_mat, theta_grad_mat,
c.set_edgecolor("face") # Reduce aliasing in output.
plt.xlabel("$h$", labelpad=0)
plt.ylabel("$y$")
cbar = plt.colorbar(label="$\\theta_j^\\star$", orientation='horizontal')
cbar = plt.colorbar(label="$\\theta^\\star$ (rad)", orientation='horizontal')
cbar.ax.set_xticklabels(cbar.ax.get_xticklabels(), rotation='45')
plt.sca(axes[1])
......@@ -97,7 +97,7 @@ def plot_peak_sweep(h_mat, y_mat, theta_star_mat, x_star_mat, theta_grad_mat,
# plt.gca().set_xlim(h_range)
# if y_range is not None:
# plt.gca().set_ylim(y_range)
# plt.colorbar(label="$\\frac{\\theta_j'(0)}{\\theta_j^\\star/x^\\star}$")
# plt.colorbar(label="$\\frac{\\theta'(0)}{\\theta^\\star/x^\\star}$")
# plt.tight_layout()
#
# plt.figure()
......@@ -115,7 +115,7 @@ def plot_peak_sweep(h_mat, y_mat, theta_star_mat, x_star_mat, theta_grad_mat,
plt.yscale('log')
plt.scatter(y_mat, theta_star_mat, c=h_mat)
plt.xlabel('$y$')
plt.ylabel('$\\theta_j^\\star$')
plt.ylabel('$\\theta^\\star$ (rad)')
plt.colorbar(label="$h$")
plt.figure()
......@@ -123,7 +123,7 @@ def plot_peak_sweep(h_mat, y_mat, theta_star_mat, x_star_mat, theta_grad_mat,
plt.yscale('log')
plt.scatter(h_mat, theta_star_mat, c=y_mat)
plt.xlabel('$h$')
plt.ylabel('$\\theta_j^\\star$')
plt.ylabel('$\\theta^\\star$ (rad)')
plt.colorbar(label="$y$")
plt.figure()
......@@ -148,7 +148,7 @@ def plot_peak_sweep(h_mat, y_mat, theta_star_mat, x_star_mat, theta_grad_mat,
if __name__ == "__main__":
n = 20000
n_points = 16
save_file = open(f"model_outputs/peak_sweep_{n}_{n_points}x{n_points}_plate_100.csv", 'r')
save_file = open(f"../model_outputs/peak_sweep/peak_sweep_{n}_{n_points}x{n_points}_plate_100.csv", 'r')
all_hs = []
all_ys = []
......
......@@ -2,7 +2,7 @@ import numpy as np
import scipy
import os
from numerical.models.slot_opt import find_slot_peak
from numerical.models.slot.slot_opt import find_slot_peak
import numerical.util.gen_utils as gen
import numerical.bem as bem
......@@ -21,7 +21,8 @@ n = 20000
w_thresh = 15
density_ratio = 0.15
save_file = open(f"model_outputs/peak_sweep_{n}_{n_points}x{n_points}_plate_500.csv", 'a')
os.makedirs("../model_outputs/peak_sweep/", exist_ok=True)
save_file = open(f"../model_outputs/peak_sweep/peak_sweep_{n}_{n_points}x{n_points}_plate_500.csv", 'a')
for H in Hs:
if last_H is not None and (H < last_H or (np.isclose(H, last_H) and np.isclose(last_Y, np.max(Ys)))):
......@@ -37,15 +38,15 @@ for H in Hs:
if last_H is not None and last_Y is not None and np.isclose(H, last_H) and \
(Y < last_Y or np.isclose(Y, last_Y)):
continue # Skip already completed positions in the last completed slot.
_, X, theta_j, _ = find_slot_peak(W, Y, H, actual_n,
varied_slot_density_ratio=density_ratio, density_w_thresh=w_thresh,
centroids=centroids, normals=normals, areas=areas, R_inv=R_inv)
_, X, theta, _ = find_slot_peak(W, Y, H, actual_n,
varied_slot_density_ratio=density_ratio, density_w_thresh=w_thresh,
centroids=centroids, normals=normals, areas=areas, R_inv=R_inv)
center_x = 0.05
near_center_vel, _ = bem.get_jet_dir_and_sigma([center_x, Y, 0], centroids, normals, areas, R_inv=R_inv)
center_grad = np.arctan2(near_center_vel[1], near_center_vel[0]) + np.pi / 2
save_file.write(f"{H / W},{Y / W},{theta_j},{2 * X / W},{center_grad / center_x}\n")
save_file.write(f"{H / W},{Y / W},{theta},{2 * X / W},{center_grad / center_x}\n")
# Make sure the data is saved to disk every time.
save_file.flush()
......
......@@ -10,16 +10,17 @@ import common.util.file_utils as file
import scipy.sparse
# Xs = [8]
W = 4.2
H = 11.47
Ys = [2.43, 3.43]
W = 2
H = 2
Ys = [2]
x_lim = 5
Xs = np.linspace(-x_lim * W / 2, x_lim * W / 2, 100)
x_lim = 20
Xs = np.linspace(-x_lim * W / 2, x_lim * W / 2, 128)
save_to_file = True
save_to_file = False
calculate_error = False
show_plot = False
normalise = False
show_plot = True
if not (save_to_file or show_plot):
raise ValueError("No output selected for model.")
......@@ -27,15 +28,15 @@ if not (save_to_file or show_plot):
m_b = 1
n = 20000
density_ratio = 0.25
w_thresh = x_lim
w_thresh = 1.2 * x_lim
out_dir = "model_outputs/exp_comparisons"
out_dir = "../model_outputs/slot/"
if not os.path.exists(out_dir) and save_to_file:
os.makedirs(out_dir)
# centroids, normals, areas = gen.gen_slot(n=n, h=h, w=w, length=50, depth=50)
centroids, normals, areas = gen.gen_varied_slot(n=n, H=H, W=W, length=50, depth=50, w_thresh=w_thresh,
centroids, normals, areas = gen.gen_varied_slot(n=n, H=H, W=W, length=100, depth=50, w_thresh=w_thresh,
density_ratio=density_ratio)
print("Requested n = {0}, using n = {1}.".format(n, len(centroids)))
# plot_3d_point_sets([centroids])
......@@ -48,15 +49,15 @@ print(f"Condition numbers: 1 norm = {condition_number_1}, inf norm = {condition_
for Y in Ys:
if calculate_error:
theta_js = []
thetas = []
for X in Xs:
print("Testing X =", X)
R_b = bem.get_R_vector([X, Y, 0], centroids, normals)
res_vel, sigma = bem.get_jet_dir_and_sigma([X, Y, 0], centroids, normals, areas, m_0=m_b, R_inv=R_inv,
R_b=R_b)
theta_j = math.atan2(res_vel[1], res_vel[0]) + math.pi / 2
theta_js.append(theta_j)
print(" theta_j =", theta_j)
theta = math.atan2(res_vel[1], res_vel[0]) + math.pi / 2
thetas.append(theta)
print(" theta =", theta)
residual = np.abs(m_b * R_b + np.dot(R_matrix, sigma))
max_err = condition_number_inf * np.linalg.norm(residual, np.inf) / np.linalg.norm(m_b * R_b, np.inf)
......@@ -68,30 +69,36 @@ for Y in Ys:
points[:, 1] = Y
points[:, 2] = 0
vels = bem.get_jet_dirs(points, centroids, normals, areas, m_b, R_inv, verbose=True)
theta_js = np.arctan2(vels[:, 1], vels[:, 0]) + 0.5 * np.pi
thetas = np.arctan2(vels[:, 1], vels[:, 0]) + 0.5 * np.pi
xs = Xs / (0.5 * W)
if normalise:
theta_star, x_star = sorted(zip(thetas, xs))[-1]
xs = xs / x_star
thetas = thetas / theta_star
if show_plot:
fig = plt.figure()
fig.patch.set_facecolor('white')
ax = plt.gca()
ax.plot(xs, theta_js)
ax.plot(xs, thetas)
ax.set_xlabel("$x$")
ax.set_ylabel("$\\theta_j$")
ax.set_ylabel("$\\theta$ (rad)")
ax.axvline(x=-1, linestyle='--', color='gray')
ax.axvline(x=1, linestyle='--', color='gray')
plt.show()
if save_to_file:
arr = []
for i in range(len(xs)):
arr.append([xs[i], theta_js[i]])
file_path = f"{out_dir}/W{W:.2f}H{H:.2f}Y{Y:.2f}_bem_slot_prediction_{n}_{density_ratio}_{w_thresh}.csv"
file_path = f"W{W:.2f}H{H:.2f}Y{Y:.2f}_bem_slot_prediction_{n}_{density_ratio}_{w_thresh}.csv" \
if not normalise \
else f"W{W:.2f}H{H:.2f}Y{Y:.2f}_bem_slot_prediction_{n}_{density_ratio}_{w_thresh}_normalised.csv"
alph = "abcdefgh"
i = 0
while os.path.exists(file_path):
file_path = f"{out_dir}/W{W:.2f}H{H:.2f}Y{Y:.2f}_bem_slot_prediction_{n}_{density_ratio}_{w_thresh}" \
while os.path.exists(out_dir + file_path):
file_path = f"W{W:.2f}H{H:.2f}Y{Y:.2f}_bem_slot_prediction_{n}_{density_ratio}_{w_thresh}" \
f"{alph[i]}.csv"
i += 1
file.lists_to_csv("", file_path, [xs, thetas], headers=["x", "theta"])
headers = ["x", "theta"] if not normalise else ["x_hat", "theta_hat"]
file.lists_to_csv(out_dir, file_path, [xs, thetas], headers=headers)
......@@ -17,7 +17,7 @@ Y = 2
normalize = True
m_0 = 1
n = 2000
n = 20000
xs = Xs / (0.5 * W)
......@@ -36,7 +36,7 @@ for i, H in enumerate(Hs):
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
f_dir = "model_outputs/slot_h_collapse/"
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"])
......@@ -20,7 +20,7 @@ Y = 2
normalize = True
m_0 = 1
n = 2000
n = 20000
fig_width = 5.31445
fig = plt.figure(figsize=(fig_width, fig_width / 2))
......@@ -40,7 +40,7 @@ xs = Xs / (0.5 * W)
for i, H in enumerate(Hs):
print(f"Reading H = {H}")
f_dir = "model_outputs/slot_h_collapse/"
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)
......@@ -86,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/h_sweep_plot.pdf')
plt.savefig('../model_outputs/h_sweep_plot.pdf')
plt.show()
......@@ -2,7 +2,7 @@ import numpy as np
import matplotlib.pyplot as plt
from common.util.plotting_utils import initialize_plt
from numerical.models.slot_opt import find_slot_peak
from numerical.models.slot.slot_opt import find_slot_peak
Hs = np.linspace(1, 10, 10)
W = 2
......@@ -23,9 +23,9 @@ fig = plt.figure()
ax = fig.gca()
ax.plot(Hs / W, theta_stars, "k")
ax.set_xlabel("$h$")
ax.set_ylabel("$\\theta_j^\\star$")
ax.set_ylabel("$\\theta_j^\\star$ (rad)")
fig.tight_layout()
fig.savefig('model_outputs/theta_j_star_h_var.svg')
fig.savefig('../model_outputs/theta_j_star_h_var.svg')
fig = plt.figure()
ax = fig.gca()
......@@ -34,6 +34,6 @@ ax.plot(Hs / W, x_stars, "k")
ax.set_xlabel("$h$")
ax.set_ylabel("$x^\\star$")
fig.tight_layout()
fig.savefig('model_outputs/x_star_h_var.svg')
fig.savefig('../model_outputs/x_star_h_var.svg')
plt.show()
......@@ -19,7 +19,7 @@ us = []
vs = []
speeds = []
file = open(f"model_outputs/slot_vel_data/vel_sweep_n{n}_W{W:.2f}_H{H:.2f}"
file = open(f"../model_outputs/slot_vel_data/vel_sweep_n{n}_W{W:.2f}_H{H:.2f}"
f"_drat{density_ratio}_wthresh{w_thresh}_len{length}_N{N}.csv", 'r')
for line in file.readlines():
......@@ -30,6 +30,8 @@ for line in file.readlines():
vs.append(float(v))
speeds.append(np.linalg.norm([float(u), float(v)]))
Xs, Ys, us, vs, speeds = zip(*sorted(zip(Xs, Ys, us, vs, speeds), key=lambda q: (q[1], q[0])))
if len(Xs) == N ** 2:
X_mat = np.reshape(np.array(Xs), (N, N))
Y_mat = np.reshape(np.array(Ys), (N, N))
......@@ -59,7 +61,7 @@ else:
c_xs = []
c_ys = []
c_zs = []
centroids_file = open(f"model_outputs/slot_vel_data/centroids_n{n}_W{W:.2f}_H{H:.2f}"
centroids_file = open(f"../model_outputs/slot_vel_data/centroids_n{n}_W{W:.2f}_H{H:.2f}"
f"_drat{density_ratio}_wthresh{w_thresh}_len{length}.csv", 'r')
for line in centroids_file.readlines():
split = line.split(",")
......@@ -70,7 +72,7 @@ for line in centroids_file.readlines():
n_xs = []
n_ys = []
n_zs = []
normals_file = open(f"model_outputs/slot_vel_data/normals_n{n}_W{W:.2f}_H{H:.2f}"
normals_file = open(f"../model_outputs/slot_vel_data/normals_n{n}_W{W:.2f}_H{H:.2f}"
f"_drat{density_ratio}_wthresh{w_thresh}_len{length}.csv", 'r')
for line in normals_file.readlines():
split = line.split(",")
......@@ -92,9 +94,9 @@ plot_log = False
if plot_log:
fig = plt.figure()
fig.gca().set_aspect('equal', 'box')
min_q_idx = 0
cnt = plt.contourf((2 * X_mat / W)[min_q_idx:, :], (Y_mat / W)[min_q_idx:, :],
(np.log(np.abs(np.arctan2(V, U) + np.pi / 2)))[min_q_idx:, :], levels=128,
min_y_idx = 0
cnt = plt.contourf((2 * X_mat / W)[min_y_idx:, :], (Y_mat / W)[min_y_idx:, :],
(np.log(np.abs(np.arctan2(V, U) + np.pi / 2)))[min_y_idx:, :], levels=128,
cmap=plt.cm.get_cmap('seismic'))
# plt.scatter((2 * P / w)[min_q_idx:, :], (Q / w)[min_q_idx:, :], color='k', marker='.')
plt.xlabel("$x$")
......@@ -110,7 +112,7 @@ if plot_log:
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes('right', size='3%', pad=0.1)
plt.colorbar(cnt, label="$log(|\\theta_j|)$", cax=cax)
plt.colorbar(cnt, label="$log(|\\theta|)$", cax=cax)
plt.tight_layout()
###############
......@@ -118,20 +120,23 @@ if plot_log:
###############
# NOTE: In the main plot (ax) all y values are doubled and then corrected in the tick labels to make the scaling work
# correctly. Other methods cause various axes to become the wrong sizes.
fig, ax = plt.subplots(figsize=(5.31445, 5.2))
fig, ax = plt.subplots(figsize=(5.31445, 4.3))
ax.set_aspect(1, 'box')
min_q_idx = 2
cnt = plt.contourf((2 * X_mat / W)[min_q_idx:, :], (2 * Y_mat / W)[min_q_idx:, :],
(np.arctan2(V, U) + np.pi / 2)[min_q_idx:, :], levels=64,
min_y_idx = 2
max_y_idx = -8
cnt = plt.contourf((2 * X_mat / W)[min_y_idx:max_y_idx, :], (2 * Y_mat / W)[min_y_idx:max_y_idx, :],
(np.arctan2(V, U) + np.pi / 2)[min_y_idx:max_y_idx, :], levels=64,
cmap=plt.cm.get_cmap('seismic'))
for c in cnt.collections:
c.set_edgecolor("face") # Reduce aliasing in output.
# plt.scatter((2 * P / w)[min_q_idx:, :], (Q / w)[min_q_idx:, :], color='k', marker='.')
plt.ylabel("$y$")
plt.xlim((2 * min(Xs) / W, 2 * max(Xs) / W))
plt.ylim(-2 * H / W - 0.1, 2 * max(Ys) / W)
plt.plot([min(Xs) * 2 / W, -W / 2, -W / 2, W / 2, W / 2, max(Xs) * 2 / W], [0, 0, -2 * H / W, -2 * H / W, 0, 0], 'k')
ax.set_yticklabels(ax.get_yticks() / 2) # Hack to get scaling and sizing right.
plt.ylim(-2 * H / W - 0.1, 2 * np.max(Y_mat[min_y_idx:max_y_idx, :]) / W)
plt.plot([min(Xs) * 2 / W, -1, -1, 1, 1, max(Xs) * 2 / W], [0, 0, -2 * H / W, -2 * H / W, 0, 0], 'k')
Yticks = np.array(range(-2, int(np.max(Y_mat[min_y_idx:max_y_idx, :])) + 1, 2))
ax.set_yticks(Yticks)
ax.set_yticklabels(Yticks / 2) # Hack to get scaling and sizing right.
plt.scatter(np.divide(c_xs, 0.5 * W), np.divide(c_ys, 0.5 * W), marker='.', color='k')
# for x, y, angle in zip(xs, ys, angles):
# plt.scatter(x, y, marker=(3, 0, np.degrees(angle)), color='k')
......@@ -147,15 +152,18 @@ plt.annotate(f"($a$)", xy=(0, 0), xytext=(0.025, 0.975),
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='3%', pad=0.1)
plt.colorbar(cnt, label="$\\theta_j$ (rad)", cax=cax, ticks=[-np.pi / 4, -np.pi / 8, 0, np.pi / 8, np.pi / 4])
# plt.colorbar(cnt, label="$\\theta$ (rad)", cax=cax, ticks=[-np.pi / 4, -np.pi / 8, 0, np.pi / 8, np.pi / 4])
cbar = plt.colorbar(cnt, label="$\\theta$ (rad)", cax=cax, ticks=[-np.pi / 4, -np.pi / 8, 0, np.pi / 8, np.pi / 4])
cax.set_yticklabels(["$-\\pi / 4$", "$-\\pi / 8$", "0", "$\\pi / 8$", "$\\pi / 4$"])
# cbar = plt.colorbar(cnt, label="$\\theta$ (rad)", cax=cax, ticks=[-np.pi / 8, 0, np.pi / 8])
# cax.set_yticklabels(["$-\\pi / 8$", "0", "$\\pi / 8$"])
lax = divider.append_axes('bottom', 1.2, pad=0.1, sharex=ax)
theta_js = np.arctan2(vs, us) + np.pi / 2
line_Xs, line_Ys, line_theta_js = zip(*[(X, Y, theta_j) for X, Y, theta_j in zip(Xs, Ys, theta_js) if Y / W == line_y])
lax.plot(2 * np.divide(line_Xs, W), line_theta_js, 'k--', label=f"$y = {line_y:.2f}$")
thetas = np.arctan2(vs, us) + np.pi / 2
line_Xs, line_Ys, line_thetas = zip(*[(X, Y, theta) for X, Y, theta in zip(Xs, Ys, thetas) if Y / W == line_y])
lax.plot(2 * np.divide(line_Xs, W), line_thetas, 'k--', label=f"$y = {line_y:.2f}$")
plt.xlabel("$x$")
plt.ylabel("$\\theta_j$")
plt.ylabel("$\\theta$ (rad)")
plt.legend(frameon=False, loc='lower right')
plt.annotate(f"($b$)", xy=(0, 0), xytext=(0.025, 0.95),
textcoords='axes fraction',
......
......@@ -84,7 +84,7 @@ if __name__ == "__main__":
ax2 = fig.add_subplot(312)
ax2.plot(ns, theta_js, 'k')
ax2.set_ylabel("$\\theta_j^\\star$")
ax2.set_ylabel("$\\theta_j^\\star$ (rad)")
ax2.set_xlabel("$N$")
ax2_conv = ax2.twinx()
ax2_conv.plot(ns, theta_js_conv, 'k--')
......
......@@ -8,45 +8,44 @@ import scipy.sparse
import numerical.bem as bem
import numerical.util.gen_utils as gen
if not os.path.exists("model_outputs/slot_vel_data"):
os.makedirs("model_outputs/slot_vel_data")
if not os.path.exists("../model_outputs/slot_vel_data"):
os.makedirs("../model_outputs/slot_vel_data")
offset = 0.05
W = 2
H = 2
W = 10
H = 1
N = 64
Ys = np.concatenate([np.linspace(offset, 3, np.round(3 * N / 4) - 1), [1], np.linspace(3 + 0.1, W * 5, np.ceil(N / 4))])
Ys = np.concatenate([np.linspace(offset, 3, np.round(3 * N / 4) - 1), [2], np.linspace(3 + 0.1, W * 5, np.ceil(N / 4))])
Ys = sorted(Ys)
Xs = np.concatenate([np.linspace(-7 * W / 2, -W, np.ceil(N / 8)),
x_limit = 1.5
Xs = np.concatenate([np.linspace(-x_limit * W / 2, -W, np.ceil(N / 8)),
np.linspace(-W + 0.1, W - 0.1, np.round(3 * N / 4)),
np.linspace(W, 7 * W / 2, np.ceil(N / 8))])
# qs = np.linspace(3 * w, 5 * w, 16)
# ps = np.linspace(0, 5 * w / 2, 16)
np.linspace(W, x_limit * W / 2, np.ceil(N / 8))])
Xs = sorted(Xs)
m_0 = 1
n = 20000
density_ratio = 0.1
w_thresh = 12
w_thresh = x_limit * 2
length = 100
# centroids, normals, areas = gen.gen_slot(n=n, h=h, w=w, length=50, depth=50)
centroids, normals, areas = gen.gen_varied_slot(n=n, H=H, W=W, length=length, depth=50, w_thresh=w_thresh,
density_ratio=density_ratio)
centroids_file = open(f"model_outputs/slot_vel_data/centroids_n{n}_W{W:.2f}_H{H:.2f}"
centroids_file = open(f"../model_outputs/slot_vel_data/centroids_n{n}_W{W:.2f}_H{H:.2f}"
f"_drat{density_ratio}_wthresh{w_thresh}_len{length}.csv", 'w')
for c in centroids:
centroids_file.write(f"{c[0]},{c[1]},{c[2]}\n")
centroids_file.close()
normals_file = open(f"model_outputs/slot_vel_data/normals_n{n}_W{W:.2f}_H{H:.2f}"
normals_file = open(f"../model_outputs/slot_vel_data/normals_n{n}_W{W:.2f}_H{H:.2f}"
f"_drat{density_ratio}_wthresh{w_thresh}_len{length}.csv", 'w')
for normal in normals:
normals_file.write(f"{normal[0]},{normal[1]},{normal[2]}\n")
normals_file.close()
output_path = f"model_outputs/slot_vel_data/vel_sweep_n{n}_W{W:.2f}_H{H:.2f}" \
output_path = f"../model_outputs/slot_vel_data/vel_sweep_n{n}_W{W:.2f}_H{H:.2f}" \
f"_drat{density_ratio}_wthresh{w_thresh}_len{length}_N{N}.csv"
if os.path.exists(output_path):
print("Output path already exists!")
......@@ -62,7 +61,7 @@ us = []
vs = []
for Y, X in itertools.product(Ys, Xs):
print(f"Testing p={X:5.3f}, q={Y:5.3f}")
print(f"Testing X={X:5.3f}, Y={Y:5.3f}")
R_b = bem.get_R_vector([X, Y, 0], centroids, normals)
res_vel, sigma = bem.get_jet_dir_and_sigma([X, Y, 0], centroids, normals, areas, m_0=m_0, R_inv=R_inv, R_b=R_b)
speed = np.linalg.norm(res_vel)
......@@ -79,12 +78,3 @@ P, Q = np.meshgrid(Xs, Ys)
S = np.reshape(np.array(speeds), P.shape)
U = np.reshape(np.array(us), P.shape)
V = np.reshape(np.array(vs), P.shape)
fig = plt.figure()
fig.gca().set_aspect('equal', 'box')
cnt = plt.contourf(2 * P / W, Q / W, np.abs(np.arctan2(V, U) + np.pi / 2), levels=128)
plt.xlabel("$\\bar{p}$")
plt.ylabel("$q / w$")
plt.colorbar(label="$|\\theta_j|$")
plt.plot([min(Xs), -W / 2, -W / 2, W / 2, W / 2, max(Xs)], [0, 0, -H, -H, 0, 0], 'k')
plt.show()
......@@ -37,7 +37,7 @@ for i, Y in enumerate(Ys):
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_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"])
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
from common.util.plotting_utils import initialize_plt
initialize_plt()
......@@ -18,7 +12,7 @@ Xs = np.linspace(-3 * W, 3 * W, 300)
Ys = np.linspace(1, 5, 5)
m_0 = 1
n = 2000
n = 20000
fig_width = 5.31445
fig = plt.figure(figsize=(fig_width, fig_width / 2))
......@@ -38,7 +32,7 @@ xs = Xs / (0.5 * W)
for i, Y in enumerate(Ys):
print(f"Reading Y = {Y}")
f_dir = "model_outputs/slot_y_collapse/"
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)
......