Commit d044286c authored by Elijah Andrews's avatar Elijah Andrews

Added step models.

parent 7ae50866
import numpy as np
import math
from scipy.optimize import curve_fit
import numerical.util.gen_utils as gen
import numerical.bem as bem
import matplotlib.pyplot as plt
from common.util.plotting_utils import plot_3d_point_sets
import os
import common.util.file_utils as file
import scipy.sparse
# Xs = [8]
H = 2
Ys = [1, 2, 4, 8]
x_lim = 15
Xs = np.linspace(-x_lim, x_lim, 200)
calculate_error = False
m_b = 1
n = 15000
density_ratio = 0.25
thresh_dist = x_lim
out_dir = "model_outputs/exp_comparisons"
centroids, normals, areas = gen.gen_varied_step(n=n, H=H, length=100, depth=50, thresh_dist=thresh_dist * 1.5,
density_ratio=density_ratio)
print("Requested n = {0}, using n = {1}.".format(n, len(centroids)))
plot_3d_point_sets([centroids])
R_matrix = bem.get_R_matrix(centroids, normals, areas, dtype=np.float32)
R_inv = scipy.linalg.inv(R_matrix)
condition_number_1 = np.linalg.norm(R_inv, 1) * np.linalg.norm(R_matrix, 1)
condition_number_inf = np.linalg.norm(R_inv, np.inf) * np.linalg.norm(R_matrix, np.inf)
print(f"Condition numbers: 1 norm = {condition_number_1}, inf norm = {condition_number_inf}")
fig = plt.figure()
fig.patch.set_facecolor('white')
ax = plt.gca()
for Y in Ys:
if calculate_error:
theta_js = []
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)
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)
print(f" Max res = {np.max(residual):.3e}, Mean res = {np.mean(residual):.3e},"
f" Max err = {max_err:.3e}")
else:
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_b, R_inv, verbose=True)
theta_js = np.arctan2(vels[:, 1], vels[:, 0]) + 0.5 * np.pi
ax.plot(Xs, theta_js, label=f"$Y={Y}$")
ax.set_xlabel("$X$")
ax.set_ylabel("$\\theta_j$")
ax.axvline(x=0, linestyle='--', color='gray')
ax.legend()
plt.show()
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
import numpy as np
from common.util.plotting_utils import initialize_plt
n = 15000
H = 2
density_ratio = 0.1
thresh_dist = 12
length = 100
N = 32
Xs = []
Ys = []
us = []
vs = []
speeds = []
file = open(f"model_outputs/step_vel_data/vel_sweep_n{n}_H{H:.2f}"
f"_drat{density_ratio}_thresh{thresh_dist}_len{length}_N{N}.csv", 'r')
for line in file.readlines():
X, Y, u, v = line.split(',')
Xs.append(float(X))
Ys.append(float(Y))
us.append(float(u))
vs.append(float(v))
speeds.append(np.linalg.norm([float(u), float(v)]))
if len(Xs) == N ** 2:
X_mat = np.reshape(np.array(Xs), (N, N))
Y_mat = np.reshape(np.array(Ys), (N, N))
S = np.reshape(np.array(speeds), (N, N))
U = np.reshape(np.array(us), X_mat.shape)
V = np.reshape(np.array(vs), X_mat.shape)
else:
X_mat = np.empty((N, N))
X_mat.fill(np.nan)
Y_mat = np.empty((N, N))
Y_mat.fill(np.nan)
U = np.empty((N, N))
U.fill(np.nan)
V = np.empty((N, N))
V.fill(np.nan)
S = np.empty((N, N))
S.fill(np.nan)
for k in range(len(Xs)):
i = int(np.floor(k / N))
j = int(k % N)
X_mat[i, j] = Xs[k]
Y_mat[i, j] = Ys[k]
U[i, j] = us[k]
V[i, j] = vs[k]
S[i, j] = speeds[k]
c_xs = []
c_ys = []
c_zs = []
centroids_file = open(f"model_outputs/step_vel_data/centroids_n{n}_H{H:.2f}"
f"_drat{density_ratio}_thresh{thresh_dist}_len{length}.csv", 'r')
for line in centroids_file.readlines():
split = line.split(",")
c_xs.append(float(split[0]))
c_ys.append(float(split[1]))
c_zs.append(float(split[2]))
n_xs = []
n_ys = []
n_zs = []
normals_file = open(f"model_outputs/step_vel_data/normals_n{n}_H{H:.2f}"
f"_drat{density_ratio}_thresh{thresh_dist}_len{length}.csv", 'r')
for line in normals_file.readlines():
split = line.split(",")
n_xs.append(float(split[0]))
n_ys.append(float(split[1]))
n_zs.append(float(split[2]))
min_z = np.min([z for z, x in zip(c_zs, c_xs) if min(Xs) <= x <= max(Xs)])
c_xs, c_ys, c_zs, n_xs, n_ys, n_zs = zip(*[o for o in zip(c_xs, c_ys, c_zs, n_xs, n_ys, n_zs)
if o[2] == min_z and min(Xs) <= o[0] <= max(Xs)])
angles = np.arctan2(n_ys, n_xs) - np.pi / 2
initialize_plt()
############
# LOG PLOT #
############
plot_log = False
if plot_log:
fig = plt.figure()
fig.gca().set_aspect('equal', 'box')
min_q_idx = 0
cnt = plt.contourf(X_mat[min_q_idx:, :], Y_mat[min_q_idx:, :],
(np.log(np.abs(np.arctan2(V, U) + np.pi / 2)))[min_q_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$")
plt.ylabel("$y$")
plt.xlim((min(Xs), max(Xs)))
plt.ylim(-H - 0.1, max(Ys))
plt.plot([min(Xs), 0, 0, max(Xs)], [0, 0, -H, -H], 'k')
plt.scatter(c_xs, c_ys, marker='.', color='k')
# for x, y, angle in zip(xs, ys, angles):
# plt.scatter(x, y, marker=(3, 0, np.degrees(angle)), color='k')
# plt.axhline(0.5, color='k', linestyle='--')
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.tight_layout()
###############
# NORMAL PLOT #
###############
# 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))
ax.set_aspect(1, 'box')
min_q_idx = 1
min_x_idx = 1
max_x_idx = -1
cnt = plt.contourf(X_mat[min_q_idx:, min_x_idx:max_x_idx], Y_mat[min_q_idx:, min_x_idx:max_x_idx],
(np.arctan2(V, U) + np.pi / 2)[min_q_idx:, min_x_idx:max_x_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((min(Xs), max(Xs)))
plt.xlim((-10, 10)) # TODO: Don't do this
plt.ylim(-H - 0.1, max(Ys))
plt.plot([min(Xs), 0, 0, max(Xs)], [0, 0, -H, -H], 'k')
plt.scatter(c_xs, c_ys, marker='.', color='k')
# for x, y, angle in zip(xs, ys, angles):
# plt.scatter(x, y, marker=(3, 0, np.degrees(angle)), color='k')
plt.tick_params(axis='x', which='both', labelbottom=False)
target_y = 1
Y_idx = int(np.argmin(np.abs(np.subtract(Ys, target_y))))
line_y = Ys[Y_idx]
plt.axhline(line_y, color='k', linestyle='--')
plt.annotate(f"($a$)", xy=(0, 0), xytext=(0.025, 0.975),
textcoords='axes fraction',
horizontalalignment='left', verticalalignment='top')
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])
cax.set_yticklabels(["$-\\pi / 4$", "$-\\pi / 8$", "0", "$\\pi / 8$", "$\\pi / 4$"])
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 == line_y])
lax.plot(line_Xs, line_theta_js, 'k--', label=f"$Y = {line_y:.2f}$")
plt.xlabel("$X$")
plt.ylabel("$\\theta_j$")
plt.legend(frameon=False, loc='lower right')
plt.annotate(f"($b$)", xy=(0, 0), xytext=(0.025, 0.95),
textcoords='axes fraction',
horizontalalignment='left', verticalalignment='top')
plt.tight_layout()
plt.show()
import numpy as np
import math
from scipy.optimize import curve_fit
import numerical.util.gen_utils as gen
import numerical.bem as bem
import matplotlib.pyplot as plt
import itertools
from common.util.plotting_utils import plot_3d_point_sets
import os
import common.util.file_utils as file
import scipy.sparse
offset = 0.2
W = 2
H = 4
N = 16
Ys = np.concatenate([np.linspace(offset, 2, np.round(3 * N / 4)), np.linspace(2 + 0.1, 4, np.ceil(N / 4))])
Xs = np.concatenate([np.linspace(-3 * W / 2, -W, np.ceil(N / 8)),
np.linspace(-W + 0.1, W - 0.1, np.round(3 * N / 4)),
np.linspace(W, 3 * W / 2, np.ceil(N / 8))])
inner_Ys = np.linspace(offset, -H + offset, np.ceil(N / 2))
inner_Xs = np.linspace(-W / 2 + offset, W / 2 - offset, np.ceil(N / 2))
m_0 = 1
n = 15000
density_ratio = 0.25
w_thresh = 6
# 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,
density_ratio=density_ratio)
print("Requested n = {0}, using n = {1}.".format(n, len(centroids)))
# plot_3d_point_sets([centroids])
R_matrix = bem.get_R_matrix(centroids, normals, areas, dtype=np.float32)
R_inv = scipy.linalg.inv(R_matrix)
speeds = []
us = []
vs = []
for Y, X in itertools.product(Ys, Xs):
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)
speeds.append(speed)
us.append(res_vel[0] / speed)
vs.append(res_vel[1] / speed)
X_mat, Y_mat = np.meshgrid(Xs, Ys)
S = np.reshape(np.array(speeds), X_mat.shape)
U = np.reshape(np.array(us), X_mat.shape)
V = np.reshape(np.array(vs), X_mat.shape)
inner_speeds = []
inner_us = []
inner_vs = []
for Y, X in itertools.product(inner_Ys, inner_Xs):
print(f"Testing p={X:5.3f}, q={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)
inner_speeds.append(speed)
inner_us.append(res_vel[0] / speed)
inner_vs.append(res_vel[1] / speed)
inner_X, inner_Y = np.meshgrid(inner_Xs, inner_Ys)
inner_S = np.reshape(np.array(inner_speeds), inner_X.shape)
inner_U = np.reshape(np.array(inner_us), inner_X.shape)
inner_V = np.reshape(np.array(inner_vs), inner_X.shape)
fig = plt.figure()
fig.gca().set_aspect('equal', 'box')
cnt = plt.contourf(X_mat, Y_mat, S, levels=32)
plt.contourf(inner_X, inner_Y, inner_S, levels=32, vmin=cnt.get_clim()[0], vmax=cnt.get_clim()[1])
plt.quiver(X_mat, Y_mat, U, V)
# plt.quiver(inner_X, inner_Y, inner_U, inner_V)
plt.xlabel("$X$")
plt.ylabel("$Y$")
plt.colorbar(label="$|v|$")
plt.plot([min(Xs), -W / 2, -W / 2, W / 2, W / 2, max(Xs)], [0, 0, -H, -H, 0, 0], 'k')
plt.show()
import itertools
import os
import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse
import numerical.bem as bem
import numerical.util.gen_utils as gen
if not os.path.exists("model_outputs/step_vel_data"):
os.makedirs("model_outputs/step_vel_data")
offset = 0.05
H = 2
N = 32
Ys = np.concatenate([np.linspace(offset, 3, np.round(3 * N / 4) - 1), [1], np.linspace(3 + 0.1, 5, np.ceil(N / 4))])
Ys = sorted(Ys)
Xs = np.concatenate([np.linspace(-15, -5.1, np.ceil(N / 8)),
np.linspace(-5, 5, np.round(3 * N / 4)),
np.linspace(5.1, 15, np.ceil(N / 8))])
# qs = np.linspace(3 * w, 5 * w, 16)
# ps = np.linspace(0, 5 * w / 2, 16)
m_0 = 1
n = 15000
density_ratio = 0.1
thresh_dist = 12
length = 100
centroids, normals, areas = gen.gen_varied_step(n=n, H=H, length=length, depth=50, thresh_dist=thresh_dist,
density_ratio=density_ratio)
centroids_file = open(f"model_outputs/step_vel_data/centroids_n{n}_H{H:.2f}"
f"_drat{density_ratio}_thresh{thresh_dist}_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/step_vel_data/normals_n{n}_H{H:.2f}"
f"_drat{density_ratio}_thresh{thresh_dist}_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/step_vel_data/vel_sweep_n{n}_H{H:.2f}" \
f"_drat{density_ratio}_thresh{thresh_dist}_len{length}_N{N}.csv"
if os.path.exists(output_path):
print("Output path already exists!")
exit()
file = open(output_path, 'w')
print("Requested n = {0}, using n = {1}.".format(n, len(centroids)))
R_matrix = bem.get_R_matrix(centroids, normals, areas, dtype=np.float32)
R_inv = scipy.linalg.inv(R_matrix)
speeds = []
us = []
vs = []
for Y, X in itertools.product(Ys, Xs):
print(f"Testing p={X:5.3f}, q={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)
speeds.append(speed)
us.append(res_vel[0] / speed)
vs.append(res_vel[1] / speed)
file.write(f"{X},{Y},{res_vel[0]},{res_vel[1]}\n")
file.flush()
os.fsync(file.fileno())
file.close()
......@@ -249,6 +249,94 @@ def gen_varied_slot(n=3000, H=3, W=2, length=50, depth=50, w_thresh=6, density_r
return centroids, normals, areas
def gen_varied_step(n=3000, H=3, length=50, depth=50, thresh_dist=6, density_ratio=0.25, save_to_files=False):
"""
Generates a step with varying panel density.
:param n: Approximate number of panels
:param H: Step height
:param length: Geometry length
:param depth: Geometry depth
:param thresh_dist: threshold at which to reduce density
:param density_ratio: Ratio of densities
:return: centroids, normals, areas
"""
if thresh_dist < 0:
raise ValueError(f"gen_varied_step w_thresh cannot be less than zero ({thresh_dist})")
if thresh_dist > length / 2:
warnings.warn("threshold distance too high.")
raise ValueError
centroids = []
normals = []
areas = []
total_surface_length = length + H
dense_surface_length = H + thresh_dist * 2
sparse_surface_length = total_surface_length - dense_surface_length
zeta = dense_surface_length / (density_ratio * sparse_surface_length)
n_dense = n * zeta / (1 + zeta)
n_sparse = n - n_dense
n_wall = int(round(n_dense * H / dense_surface_length))
n_dense_surface_boundary = int(round((n_dense - n_wall) / 2))
n_sparse_surface_boundary = int(round(n_sparse / 2))
######################
# Surface boundaries #
######################
filename = "model_outputs/step_left_sparse.csv" if save_to_files else None
p_centroids, p_normals, p_areas = gen_plane([0, -H, - depth / 2],
[0, -H, depth / 2],
[thresh_dist, -H, - depth / 2],
n_dense_surface_boundary, filename=filename)
centroids.extend(p_centroids)
normals.extend(p_normals)
areas.extend(p_areas)
filename = "model_outputs/step_left_dense.csv" if save_to_files else None
p_centroids, p_normals, p_areas = gen_plane([thresh_dist, -H, - depth / 2],
[thresh_dist, -H, depth / 2],
[length / 2, -H, - depth / 2],
n_sparse_surface_boundary, filename=filename)
centroids.extend(p_centroids)
normals.extend(p_normals)
areas.extend(p_areas)
filename = "model_outputs/step_right_dense.csv" if save_to_files else None
p_centroids, p_normals, p_areas = gen_plane([-thresh_dist, 0, - depth / 2],
[-thresh_dist, 0, depth / 2],
[0, 0, - depth / 2],
n_dense_surface_boundary, filename=filename)
centroids.extend(p_centroids)
normals.extend(p_normals)
areas.extend(p_areas)
filename = "model_outputs/step_right_sparse.csv" if save_to_files else None
p_centroids, p_normals, p_areas = gen_plane([-length / 2, 0, - depth / 2],
[-length / 2, 0, depth / 2],
[-thresh_dist, 0, - depth / 2],
n_sparse_surface_boundary, filename=filename)
centroids.extend(p_centroids)
normals.extend(p_normals)
areas.extend(p_areas)
######################
# Step wall #
######################
filename = "model_outputs/step_wall.csv" if save_to_files else None
p_centroids, p_normals, p_areas = gen_plane([0, 0, - depth / 2],
[0, 0, depth / 2],
[0, -H, - depth / 2],
n_wall, filename=filename)
centroids.extend(p_centroids)
normals.extend(p_normals)
areas.extend(p_areas)
return centroids, normals, areas
def gen_corner(n=3000, length=25, depth=50, angle=np.pi / 2):
centroids = []
normals = []
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment