diff --git a/pymead/analysis/calc_aero_data.py b/pymead/analysis/calc_aero_data.py
index e03fbcf1..f7da112a 100644
--- a/pymead/analysis/calc_aero_data.py
+++ b/pymead/analysis/calc_aero_data.py
@@ -3,6 +3,7 @@
import typing
from copy import deepcopy
import time
+from collections import namedtuple
import numpy as np
from matplotlib import pyplot as plt
@@ -11,7 +12,7 @@
from pymead.analysis.read_aero_data import read_aero_data_from_xfoil, read_Cp_from_file_xfoil, read_bl_data_from_mses, \
read_forces_from_mses, read_grid_stats_from_mses, read_field_from_mses, read_streamline_grid_from_mses, \
- flow_var_idx, convert_blade_file_to_3d_array, read_actuator_disk_data_mses
+ flow_var_idx, convert_blade_file_to_3d_array, read_actuator_disk_data_mses, read_Mach_from_mses_file
from pymead.utils.file_conversion import convert_ps_to_svg
from pymead.utils.geometry import check_airfoil_self_intersection, convert_numpy_array_to_shapely_LineString, \
calculate_area_triangle_heron
@@ -248,9 +249,10 @@ def calculate_aero_data(airfoil_coord_dir: str, airfoil_name: str, coords: typin
if mplot_settings["CPK"]:
try:
- outputs_CPK = calculate_CPK_mses(os.path.join(airfoil_coord_dir, airfoil_name))
- except:
- outputs_CPK = {"CPK": 1e9}
+ outputs_CPK = calculate_CPK_power_consumption(os.path.join(airfoil_coord_dir, airfoil_name))
+ except Exception as e:
+ print(f"{e = }")
+ outputs_CPK = {"CPK": 1e9, "diss_shock": 1e9, "diss_surf": 1e9, "Edot": 1e9}
aero_data = {**aero_data, **outputs_CPK}
t2 = time.time()
@@ -1007,10 +1009,178 @@ def line_integral_CPK_inviscid_old(Cp: np.ndarray, rho: np.ndarray, u: np.ndarra
return integral
+def line_integral_Edot_inviscid(Cp: np.ndarray, rho: np.ndarray, u: np.ndarray,
+ v: np.ndarray, V: np.ndarray, x: np.ndarray, y: np.ndarray, n_hat_dir: str):
+ """
+ ... according to
+ Drela's "Power Balance in Aerodynamic Flows", a 2009 AIAA Journal article
+
+ Parameters
+ ==========
+ Cp: np.ndarray
+ Node-centered pressure coefficient field
+
+ rho: np.ndarray
+ Node-centered density field
+
+ u: np.ndarray
+ Node-centered x-velocity field
+
+ v: np.ndarray
+ Node-centered y-velocity field
+
+ V: np.ndarray
+ Node-centered velocity magnitude field
+
+ x: np.ndarray
+ 1-d array of x-coordinates of the line along which to compute the line integral
+
+ y: np.ndarray
+ 1-d array of y-coordinates of the line along which to compute the line integral
+
+ n_hat_right: bool
+ Whether the perpendicular n_hat vector points to the right of the line (If ``True``, computed as the sum of
+ the arctan of the value of dx/dy of the line less 90 degrees. If ``False``, 90 degrees is added instead).
+ """
+ if n_hat_dir not in ["up", "down", "left", "right"]:
+ raise ValueError("Invalid direction for n_hat")
+
+ perp_angle_dict = {
+ "right": -np.pi / 2,
+ "left": np.pi / 2,
+ "up": -np.pi / 2,
+ "down": np.pi / 2
+ }
+
+ # Calculate the direction of n_hat
+ with np.errstate(divide="ignore", invalid="ignore"):
+ # if n_hat_dir in ["right", "left"]:
+ # angle = np.arctan2(np.gradient(y, x), 1)
+ # angle = np.arctan2(y[1:] - y[:-1], x[1:] - x[:-1])
+ # else:
+ angle = np.arctan2(1, np.gradient(x, y))
+
+ # temp_angle = np.zeros(angle.shape)
+ # for idx, a in enumerate(angle):
+ # if np.isnan(a):
+ # if idx == 0:
+ # temp_angle[idx] = angle[idx + 1]
+ # else:
+ # temp_angle[idx] = angle[idx - 1]
+ # else:
+ # temp_angle[idx] = angle[idx]
+ # angle = temp_angle
+
+ perp_angle = perp_angle_dict[n_hat_dir]
+ n_hat = np.column_stack((np.cos(angle + perp_angle), np.sin(angle + perp_angle)))
+ V_vec = np.column_stack((u, v))
+
+ # Compute the dot product V_vec * n_hat
+ dot_product = np.array([np.dot(V_vec_i, n_hat_i) for V_vec_i, n_hat_i in zip(V_vec, n_hat)])
+ # integrand = (-Cp - rho * V**2 + rho) * dot_product.flatten()
+ integrand = (Cp + rho * (V**2 - 1)) * dot_product.flatten()
+
+ # Build the length increment vector (dl)
+ dl = np.array([0.0])
+ dl = np.append(dl, np.hypot(x[1:] - x[:-1], y[1:] - y[:-1])) # compute incremental length along x and y
+ # print(f"Before cumulative summation, {dl = }")
+ # print(f"{x = }")
+ # print(f"{y = }")
+ dl = np.cumsum(dl)
+
+ # Integrate
+ integral = np.trapz(integrand, dl)
+ return integral
+
+
+def viscous_Edot(last_BL_point: typing.Dict[str, np.ndarray]):
+ """
+ Eq. (75) of Drela's Power Balance in Aerodynamic Flows in 2-D form
+ """
+ if any([k not in last_BL_point.keys() for k in ("rhoe/rhoinf", "Ue/Uinf", "theta", "theta*")]):
+ raise ValueError(f"Missing a boundary layer key in the Edot calculation. {last_BL_point.keys() = }")
+
+ delta_K = 2 * last_BL_point["theta"] - last_BL_point["theta*"]
+ Edot_visc = np.sum(last_BL_point["rhoe/rhoinf"] * last_BL_point["Ue/Uinf"]**3 * delta_K)
+ return Edot_visc
+
+
+def surface_dissipation(last_BL_point: typing.Dict[str, np.ndarray]):
+ """
+ Eq. (73) of Drela's Power Balance in Aerodynamic Flows
+ """
+ if any([k not in last_BL_point.keys() for k in ("K",)]):
+ raise ValueError(f"Missing a boundary layer key in the Edot calculation. {last_BL_point.keys() = }")
+
+ diss = np.sum(2 * last_BL_point["K"])
+ return diss
+
+
+def shock_dissipation(x_edge: typing.List[np.ndarray], y_edge: typing.List[np.ndarray], u: np.ndarray,
+ v: np.ndarray, dCpt: np.ndarray):
+ Shock = namedtuple("Shock", ("j", "i_start", "i_end"))
+ shock_indices = np.argwhere(dCpt < -0.002)
+ shock_indices = np.fliplr(shock_indices)
+ shock_indices = shock_indices[shock_indices[:, 1].argsort()]
+ shock_indices = shock_indices[shock_indices[:, 0].argsort(kind="mergesort")]
+ shocks = {k: [] for k in np.unique(shock_indices[:, 0])}
+
+ i_start = 0
+ for idx, row in enumerate(shock_indices):
+ j = row[0]
+ if idx == 0:
+ i_start = row[1]
+ continue
+ if row[0] == shock_indices[idx, 0] and row[1] == shock_indices[idx - 1, 1] + 1:
+ pass
+ else:
+ i_end = shock_indices[idx - 1, 1] + 1
+ shocks[j].append(Shock(j=j, i_start=i_start, i_end=i_end))
+ i_start = row[1]
+
+ i1 = 0
+ edge_split_i = []
+ for _idx, edge_array in enumerate(x_edge):
+ i2 = i1 + edge_array.shape[1] - 2
+ edge_split_i.append([i1, i2])
+ i1 = i2 + 1
+
+ def integral(s: Shock):
+ total_dCpt = np.sum(dCpt[s.i_start:s.i_end, s.j])
+ _u = u[s.i_start - 1, s.j]
+ _v = v[s.i_start - 1, s.j]
+
+ j_flow = 0
+ j_grid = 0
+ for flow_idx, split_idx in enumerate(edge_split_i):
+ if split_idx[0] <= s.j <= split_idx[1]:
+ j_flow = flow_idx
+ j_grid = s.j - split_idx[0]
+ break
+
+ x1 = x_edge[j_flow][s.i_start, j_grid]
+ x2 = x_edge[j_flow][s.i_start, j_grid + 1]
+ y1 = y_edge[j_flow][s.i_start, j_grid]
+ y2 = y_edge[j_flow][s.i_start, j_grid + 1]
+ n_hat_angle = np.arctan2(y2 - y1, x2 - x1) - np.pi / 2
+ n_hat = np.array([np.cos(n_hat_angle), np.sin(n_hat_angle)])
+ dot_product = np.dot(np.array([_u, _v]), n_hat)
+ dl = np.hypot(x2 - x1, y2 - y1)
+ shock_int = np.abs(total_dCpt) * dot_product * dl
+ return shock_int
+
+ shock_diss = 0.0
+ for shock_list in shocks.values():
+ for shock in shock_list:
+ shock_diss += integral(shock)
+
+ return shock_diss
+
+
def line_integral_CPK_inviscid(Cp_up: np.ndarray, Cp_down: np.ndarray, rho_up: np.ndarray, rho_down: np.ndarray,
- u_up: np.ndarray, u_down: np.ndarray,
- v_up: np.ndarray, v_down: np.ndarray, V_up: np.ndarray, V_down: np.ndarray,
- x: np.ndarray, y: np.ndarray):
+ u_up: np.ndarray, u_down: np.ndarray,
+ v_up: np.ndarray, v_down: np.ndarray, V_up: np.ndarray, V_down: np.ndarray,
+ x: np.ndarray, y: np.ndarray):
"""
Computes the mechanical flow power line integral in the inviscid streamtube only according to
Drela's "Power Balance in Aerodynamic Flows", a 2009 AIAA Journal article
@@ -1286,6 +1456,137 @@ def calculate_CPK_mses_old(analysis_subdir: str, configuration: str = "underwing
return {"CPK": CPK, "capSS": capSS, "epma": epma}
+def calculate_CPK_power_consumption(analysis_subdir: str):
+ """
+ A specialized function that calculates the mechanical flow power coefficient for an underwing trailing edge
+ aero-propulsive configuration.
+ """
+ airfoil_system_name = os.path.split(analysis_subdir)[-1]
+ field_file = os.path.join(analysis_subdir, f'field.{airfoil_system_name}')
+ grid_stats_file = os.path.join(analysis_subdir, 'mplot_grid_stats.log')
+ grid_file = os.path.join(analysis_subdir, f'grid.{airfoil_system_name}')
+ blade_file = os.path.join(analysis_subdir, f"blade.{airfoil_system_name}")
+ bl_file = os.path.join(analysis_subdir, f"bl.{airfoil_system_name}")
+ mses_file = os.path.join(analysis_subdir, f"mses.{airfoil_system_name}")
+ mses_log_file = os.path.join(analysis_subdir, "mses.log")
+ coords = convert_blade_file_to_3d_array(blade_file)
+
+ M_inf = read_Mach_from_mses_file(mses_file)
+ gam = 1.4 # Hard-coded specific heat ratio
+ field = read_field_from_mses(field_file, M_inf=M_inf, gam=gam)
+ bl_data = read_bl_data_from_mses(bl_file)
+ grid_stats = read_grid_stats_from_mses(grid_stats_file)
+ x_grid, y_grid = read_streamline_grid_from_mses(grid_file, grid_stats)
+
+ Edot = 0.0
+
+ start_idx, end_idx = 0, x_grid[0].shape[1] - 1
+
+ for flow_section_idx in range(grid_stats["numel"] + 1):
+ Cp = convert_cell_centered_to_edge_centered(x_grid[flow_section_idx].shape,
+ field[flow_var_idx["Cp"]][:, start_idx:end_idx])
+ rho = convert_cell_centered_to_edge_centered(x_grid[flow_section_idx].shape,
+ field[flow_var_idx["rho"]][:, start_idx:end_idx])
+ u = convert_cell_centered_to_edge_centered(x_grid[flow_section_idx].shape,
+ field[flow_var_idx["u"]][:, start_idx:end_idx])
+ v = convert_cell_centered_to_edge_centered(x_grid[flow_section_idx].shape,
+ field[flow_var_idx["v"]][:, start_idx:end_idx])
+ V = convert_cell_centered_to_edge_centered(x_grid[flow_section_idx].shape,
+ field[flow_var_idx["q"]][:, start_idx:end_idx])
+ x = x_grid[flow_section_idx]
+ y = y_grid[flow_section_idx]
+
+ # Integrate over the propulsor outlet for the given flow section
+
+ inlet_integral = line_integral_Edot_inviscid(Cp[0, :], rho[0, :], u[0, :], v[0, :], V[0, :], x[0, :],
+ y[0, :], n_hat_dir="left")
+ Edot += inlet_integral
+
+ outlet_integral = line_integral_Edot_inviscid(Cp[-1, :], rho[-1, :], u[-1, :], v[-1, :], V[-1, :], x[-1, :],
+ y[-1, :], n_hat_dir="right")
+ Edot += outlet_integral
+
+ # The side cylinder contributions are negligible
+ # if flow_section_idx == len(x_grid) - 1:
+ #
+ # top_integral = line_integral_Edot_inviscid(Cp[:, -1], rho[:, -1], u[:, -1], v[:, -1], V[:, -1], x[:, -1],
+ # y[:, -1], n_hat_dir="up")
+ # CPK += top_integral
+ # print(f"{top_integral = }")
+ #
+ # elif flow_section_idx == 0:
+ # bottom_integral = line_integral_Edot_inviscid(Cp[:, 0], rho[:, 0], u[:, 0], v[:, 0], V[:, 0], x[:, 0],
+ # y[:, 0], n_hat_dir="down")
+ # CPK += bottom_integral
+ # print(f"{bottom_integral = }")
+
+ # print(f"{outlet_integral = }, {inlet_integral = }, {Edot = }")
+
+ # print(f"After adding the inlet BL contributions, {CPK - CPK_temp = }")
+
+ if flow_section_idx < grid_stats["numel"]:
+ start_idx = end_idx
+ end_idx += x_grid[flow_section_idx + 1].shape[1] - 1
+
+ end_point_counter = 0
+ last_point_idx = -1
+ max_end_point_attempts = 10
+
+ while True:
+ if end_point_counter > max_end_point_attempts:
+ raise ValueError("Reached maximum attempts for setting the boundary layer end point.")
+
+ last_BL_point = {}
+ for side in bl_data:
+ for k in ("rhoe/rhoinf", "Ue/Uinf", "theta", "theta*", "K"):
+ if k not in last_BL_point.keys():
+ last_BL_point[k] = []
+ val = side[k][last_point_idx]
+ if val <= 1e-12 and k == "K":
+ end_point_counter += 1
+ print(f"Decrementing last_point_idx...")
+ last_point_idx -= 1
+ continue
+ last_BL_point[k].append(val)
+
+ for k, v in last_BL_point.items():
+ last_BL_point[k] = np.array(v)
+ break
+
+ print(f"{Edot = }")
+
+ Edot_visc = viscous_Edot(last_BL_point=last_BL_point)
+
+ Edot += Edot_visc
+
+ print(f"{Edot_visc = }")
+
+ surface_diss = surface_dissipation(last_BL_point=last_BL_point)
+
+ print(f"{surface_diss = }")
+
+ shock_diss = shock_dissipation(x_edge=x_grid, y_edge=y_grid, u=field[flow_var_idx["u"]][:, :],
+ v=field[flow_var_idx["v"]][:, :], dCpt=field[flow_var_idx["dCpt"]][:, :])
+
+ print(f"{shock_diss = }")
+
+ CPK = Edot + surface_diss + shock_diss
+
+ print(f"{CPK = }")
+
+ if np.isnan(CPK):
+ CPK = 1e9
+
+ # print(f"{CPK = }, {inlet_integral = }, {outlet_integral = }")
+
+ # if calculate_capSS:
+ # return CPK, capSS
+ # else:
+ # return CPK
+
+ return {"CPK": CPK, "Edot": Edot, "diss_surf": surface_diss, "diss_shock": shock_diss}
+
+
def calculate_CPK_mses(analysis_subdir: str, configuration: str = "underwing_te"):
"""
A specialized function that calculates the mechanical flow power coefficient for an underwing trailing edge
diff --git a/pymead/analysis/read_aero_data.py b/pymead/analysis/read_aero_data.py
index 86ba675c..31c2a3da 100644
--- a/pymead/analysis/read_aero_data.py
+++ b/pymead/analysis/read_aero_data.py
@@ -6,7 +6,7 @@
import numpy as np
-flow_var_idx = {'M': 7, 'Cp': 8, 'p': 3, 'rho': 2, 'u': 4, 'v': 5, 'q': 6}
+flow_var_idx = {'M': 7, 'Cp': 8, 'p': 3, 'rho': 2, 'u': 4, 'v': 5, 'q': 6, "Cpt": 9, "dCpt": 10}
def read_Cp_from_file_xfoil(fname: str):
@@ -297,7 +297,16 @@ def read_bl_data_from_mses(src_file: str):
return bl
-def read_field_from_mses(src_file: str):
+def read_Mach_from_mses_file(mses_file: str):
+ with open(mses_file, "r") as f:
+ lines = f.readlines()
+ mach_line = lines[2]
+ mach_split = mach_line.split()
+ mach_float = float(mach_split[0])
+ return mach_float
+
+
+def read_field_from_mses(src_file: str, M_inf: float = None, gam: float = None):
r"""
Reads a field dump file from MSES (by default, of the form ``field.*``) and outputs the information to an array.
The array has shape :math:`9 \times m \times n`, where :math:`m` is the number of streamlines and :math:`n`
@@ -339,6 +348,16 @@ def read_field_from_mses(src_file: str):
field_array = np.array([data[:, i].reshape(n_streamlines, n_streamwise_lines).T for i in range(n_flow_vars)])
+ if M_inf is not None and gam is not None:
+ Cpt = field_array[3][:, :] * (1 + (gam - 1) / 2 * field_array[7][:, :] ** 2) ** (gam / (gam - 1)) * (2 / gam / M_inf**2)
+ Cpt = np.array([Cpt])
+ field_array = np.concatenate((field_array, Cpt))
+ dCpt = np.zeros(Cpt[0].shape)
+ for idx in range(1, Cpt[0].shape[0]):
+ dCpt[idx, :] = Cpt[0][idx, :] - Cpt[0][idx - 1, :]
+ dCpt = np.array([dCpt])
+ field_array = np.concatenate((field_array, dCpt))
+
return field_array
diff --git a/pymead/plugins/IGES/entity.py b/pymead/plugins/IGES/entity.py
index ea19af81..ffd29fd2 100644
--- a/pymead/plugins/IGES/entity.py
+++ b/pymead/plugins/IGES/entity.py
@@ -46,6 +46,7 @@ def __init__(self, ID: int, parameter_data: typing.List[IGESParam]):
self.parameter_data = parameter_data
self.param_delimiter = None
self.record_delimiter = None
+ # TODO: make the color for each entity (or at least all entities) an option in the IGES export (QComboBox?)
def write_entity_string(self, entity_starting_line: int, data_starting_line: int, data_string_lines: int):
diff --git a/pymead/post/mses_field.py b/pymead/post/mses_field.py
index cc6dc1f6..aa8d7393 100644
--- a/pymead/post/mses_field.py
+++ b/pymead/post/mses_field.py
@@ -3,7 +3,7 @@
import numpy as np
import os
from pymead.analysis.read_aero_data import read_grid_stats_from_mses, read_field_from_mses, \
- read_streamline_grid_from_mses, read_bl_data_from_mses
+ read_streamline_grid_from_mses, read_bl_data_from_mses, read_Mach_from_mses_file
from pymead.analysis.calc_aero_data import convert_cell_centered_to_edge_centered, extrapolate_data_line_mses_field
from pymead.analysis.read_aero_data import flow_var_idx
@@ -15,7 +15,9 @@
'u': 'Velocity-x (u/V\u221e)',
'v': 'Velocity-y (v/V\u221e)',
'V': 'Velocity-mag (V/V\u221e)',
- 'q': 'Speed of Sound (q/V\u221e)'}
+ 'q': 'Speed of Sound (q/V\u221e)',
+ "Cpt": "Total Pressure / Pinf",
+ "dCpt": "Delta Total Pressure"}
def generate_field_matplotlib(axs: plt.Axes or None, analysis_subdir: str, var: str, cmap_field: mpl_colors.Colormap or str ,
@@ -25,6 +27,7 @@ def generate_field_matplotlib(axs: plt.Axes or None, analysis_subdir: str, var:
field_file = os.path.join(analysis_subdir, f'field.{os.path.split(analysis_subdir)[-1]}')
grid_stats_file = os.path.join(analysis_subdir, 'mplot_grid_stats.log')
grid_file = os.path.join(analysis_subdir, f'grid.{os.path.split(analysis_subdir)[-1]}')
+ mses_file = os.path.join(analysis_subdir, f"mses.{os.path.split(analysis_subdir)[-1]}")
# bl_file = os.path.join(analysis_subdir, f'bl.{os.path.split(analysis_subdir)[-1]}')
if not os.path.exists(field_file):
raise OSError(f"Field file {field_file} not found")
@@ -34,7 +37,9 @@ def generate_field_matplotlib(axs: plt.Axes or None, analysis_subdir: str, var:
raise OSError(f"Grid file {grid_file} not found")
# bl_data = read_bl_data_from_mses(bl_file)
- field = read_field_from_mses(field_file)
+ M_inf = read_Mach_from_mses_file(mses_file)
+ gam = 1.4
+ field = read_field_from_mses(field_file, M_inf=M_inf, gam=gam)
grid_stats = read_grid_stats_from_mses(grid_stats_file)
x_grid, y_grid = read_streamline_grid_from_mses(grid_file, grid_stats)
diff --git a/pymead/post/post_process.py b/pymead/post/post_process.py
index eab6384e..1e69ac77 100644
--- a/pymead/post/post_process.py
+++ b/pymead/post/post_process.py
@@ -6,6 +6,7 @@
import numpy
from matplotlib import animation
from matplotlib.lines import Line2D
+from matplotlib import ticker as mticker
import scienceplots
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.ticker import MultipleLocator
@@ -447,58 +448,77 @@ def generate_aero_force_plot(self, var: str or list = None):
show_save_fig(fig, save_base_dir=self.image_dir, file_name_stub=f'design_{v}')
def pareto_front(self, opt_start: int = 5, opt_end: int = None, opt_num: int = 10):
- index = self.set_index(None)
- if opt_end is None:
- opt_end = index[-1]
- all_opt_F = [self.get_full_gen_opt_F(os.path.join(self.analysis_dir, f"algorithm_gen_{i}.pkl"))
- if i != 0 else self.get_full_gen_opt_F(alg_file=None) for i in index]
- all_pop_F = [self.get_full_gen_pop_F(os.path.join(self.analysis_dir, f"algorithm_gen_{i}.pkl"))
- if i != 0 else self.get_full_gen_pop_F(alg_file=None) for i in index]
- opt_range = np.linspace(opt_start, opt_end, opt_num, endpoint=True).astype(int)
- cmap = plt.get_cmap("plasma", len(opt_range))
- cmap_dict = {v: i for i, v in enumerate(opt_range)}
- fig, axs = plt.subplots(figsize=(8, 6))
- pop_labeled = False
- opt_labeled = False
- for i, F in enumerate(all_pop_F):
- if i != 0:
- extra_kwargs = {}
- if not pop_labeled:
- extra_kwargs["label"] = "Individual"
- pop_labeled = True
- valid_rows = np.argwhere(F[:, 0] < 1.0).flatten()
- axs.plot(F[valid_rows, 0], F[valid_rows, 1], ls="none", marker="s", mfc="lightgray", mec="lightgray",
- markersize=2, **extra_kwargs)
- for i, F in enumerate(all_opt_F):
- if F.ndim > 1:
- sorted_order = np.argsort(F[:, 0], axis=None)
- F[:, 0] = F[sorted_order, 0]
- F[:, 1] = F[sorted_order, 1]
- star_kwargs = dict(ls="--", marker="*")
- gray_circle_kwargs = dict(ls="none", marker="o", mfc="none", mec="gray", markersize=3)
- base_kwargs = dict(ls="none", marker="x", mfc="red", mec="red")
- if i == 0:
- axs.plot(F[0], F[1], **base_kwargs, label="Baseline")
- else:
- if i in opt_range:
- axs.plot(F[:, 0], F[:, 1], **star_kwargs, zorder=100, label=f"Opt, Gen {i}",
- color=cmap(cmap_dict[i]),
- mfc=cmap(cmap_dict[i]),
- mec=cmap(cmap_dict[i])
- )
- else:
+
+ fig, axs = plt.subplots(figsize=(10, 6))
+
+ def generate_plots(ax: plt.Axes, opt_end_):
+ index = self.set_index(None)
+ if opt_end_ is None:
+ opt_end_ = index[-1]
+ all_opt_F = [self.get_full_gen_opt_F(os.path.join(self.analysis_dir, f"algorithm_gen_{i}.pkl"))
+ if i != 0 else self.get_full_gen_opt_F(alg_file=None) for i in index]
+ all_pop_F = [self.get_full_gen_pop_F(os.path.join(self.analysis_dir, f"algorithm_gen_{i}.pkl"))
+ if i != 0 else self.get_full_gen_pop_F(alg_file=None) for i in index]
+ opt_range = np.linspace(opt_start, opt_end_, opt_num, endpoint=True).astype(int)
+ cmap = plt.get_cmap("plasma", len(opt_range))
+ cmap_dict = {v: i for i, v in enumerate(opt_range)}
+
+ pop_labeled = False
+ opt_labeled = False
+ for i, F in enumerate(all_pop_F):
+ if i != 0:
extra_kwargs = {}
- if not opt_labeled:
- extra_kwargs["label"] = "Optimal"
- opt_labeled = True
- axs.plot(F[:, 0], F[:, 1], **gray_circle_kwargs, **extra_kwargs)
- axs.set_xlim([0.058, 0.158])
- axs.set_xlabel(r"$J_P$", fontdict=font)
- axs.set_ylabel(r"$J_F$", fontdict=font)
+ if not pop_labeled:
+ extra_kwargs["label"] = "Individual"
+ pop_labeled = True
+ valid_rows = np.argwhere(F[:, 0] < 1.0).flatten()
+ ax.plot(F[valid_rows, 0], F[valid_rows, 1], ls="none", marker="s", mfc="lightgray", mec="lightgray",
+ markersize=2, **extra_kwargs)
+ for i, F in enumerate(all_opt_F):
+ if F.ndim > 1:
+ sorted_order = np.argsort(F[:, 0], axis=None)
+ F[:, 0] = F[sorted_order, 0]
+ F[:, 1] = F[sorted_order, 1]
+ star_kwargs = dict(ls="--", marker="*")
+ gray_circle_kwargs = dict(ls="none", marker="o", mfc="none", mec="gray", markersize=3)
+ base_kwargs = dict(ls="none", marker="x", mfc="red", mec="red")
+ if i == 0:
+ axs.plot(F[0], F[1], **base_kwargs, label="Baseline")
+ else:
+ if i in opt_range:
+ ax.plot(F[:, 0], F[:, 1], **star_kwargs, zorder=100, label=f"Opt, Gen {i}",
+ color=cmap(cmap_dict[i]),
+ mfc=cmap(cmap_dict[i]),
+ mec=cmap(cmap_dict[i])
+ )
+ else:
+ extra_kwargs = {}
+ if not opt_labeled:
+ extra_kwargs["label"] = "Optimal"
+ opt_labeled = True
+ ax.plot(F[:, 0], F[:, 1], **gray_circle_kwargs, **extra_kwargs)
+ ax.set_xlim([0.058, 0.19])
+ ax.set_xlabel(r"$J_P$", fontdict=font)
+ ax.set_ylabel(r"$J_F$", fontdict=font)
+ format_axis_scientific(ax)
+
+ generate_plots(axs, opt_end)
legend_font_dict = {k: v for k, v in font.items() if k != "color"}
legend_font_dict["size"] = 12
- axs.legend(prop=legend_font_dict)
- format_axis_scientific(axs)
+ axs.legend(prop=legend_font_dict, loc="lower right", ncol=2)
+
+ axs_inset = axs.inset_axes([0.6, 0.53, 0.35, 0.4])
+ generate_plots(axs_inset, opt_end)
+ x1, x2, y1, y2 = 0.068, 0.07, 0.012, 0.014
+ axs_inset.set_xlim(x1, x2)
+ axs_inset.set_ylim(y1, y2)
+ # x_labels = [item if idx % 2 else "" for idx, item in enumerate(axs_inset.get_xticklabels())]
+ # y_labels = [item if not idx % 2 else "" for idx, item in enumerate(axs_inset.get_yticklabels())]
+ # axs_inset.set_xticklabels(x_labels)
+ # axs_inset.set_yticklabels(y_labels)
+ axs_inset.xaxis.set_major_locator(mticker.MaxNLocator(2))
+ axs_inset.yaxis.set_major_locator(mticker.MaxNLocator(2, prune="lower"))
+
show_save_fig(fig, save_base_dir=self.image_dir, file_name_stub="pareto_front")
def compare_geometries(self, index: list, plot_actuator_disk: bool = True,
@@ -790,7 +810,9 @@ def generate_single_field(self, var: str, index: int, index_list, cmap_field: mp
'rho': r'Density ($\rho/\rho_\infty$)',
'u': r'Velocity-x ($u/V_\infty$)',
'v': r'Velocity-y ($v/V_\infty$)',
- 'q': r'Speed of Sound ($q/V_\infty$)'}
+ 'q': r'Speed of Sound ($q/V_\infty$)',
+ "Cpt": r"Total Pressure Over P_inf",
+ "dCpt": r"Delta Total Pressure"}
post_process_forces = load_data(
os.path.splitext(
self.post_process_force_file)[0] + weight_str + os.path.splitext(self.post_process_force_file)[1])
@@ -901,7 +923,8 @@ def generate_field_gif_matplotlib(self, var: str,
'rho': r'Density ($\rho/\rho_\infty$)',
'u': r'Velocity-x ($u/V_\infty$)',
'v': r'Velocity-y ($v/V_\infty$)',
- 'q': r'Speed of Sound ($q/V_\infty$)'}
+ 'q': r'Speed of Sound ($q/V_\infty$)',
+ "Cpt": "Total Pressure Coefficient"}
airfoil_color = 'black'
diff --git a/pymead/tests/aero_tests/Edot_test.py b/pymead/tests/aero_tests/Edot_test.py
new file mode 100644
index 00000000..3d41c90e
--- /dev/null
+++ b/pymead/tests/aero_tests/Edot_test.py
@@ -0,0 +1,45 @@
+import unittest
+
+from pymead.analysis.calc_aero_data import calculate_CPK_power_consumption
+
+
+class EdotTest(unittest.TestCase):
+ def test_Edot_baseline(self):
+ analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_70\analysis\analysis_0_w50-50"
+ Edot = calculate_CPK_power_consumption(analysis_dir)
+ print(f"Baseline: {Edot = }")
+
+ def test_Edot_opt(self):
+ analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_70\analysis\analysis_500_w50-50"
+ Edot = calculate_CPK_power_consumption(analysis_dir)
+ print(f"Opt: {Edot = }")
+
+ def test_Edot_opt_oblique_shock_1(self):
+ analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_82\analysis\analysis_127_w0-100"
+ Edot = calculate_CPK_power_consumption(analysis_dir)
+ print(f"Opt: {Edot = }")
+
+ def test_Edot_opt_oblique_shock_2(self):
+ analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_82\analysis\analysis_127_w50-50"
+ Edot = calculate_CPK_power_consumption(analysis_dir)
+ print(f"Opt: {Edot = }")
+
+ def test_Edot_opt_oblique_shock_3(self):
+ analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_82\analysis\analysis_127_w60-40"
+ Edot = calculate_CPK_power_consumption(analysis_dir)
+ print(f"Opt: {Edot = }")
+
+ def test_Edot_opt_oblique_shock_4(self):
+ analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_82\analysis\analysis_127_w70-30"
+ Edot = calculate_CPK_power_consumption(analysis_dir)
+ print(f"Opt: {Edot = }")
+
+ def test_Edot_opt_oblique_shock_5(self):
+ analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_82\analysis\analysis_127_w100-0"
+ Edot = calculate_CPK_power_consumption(analysis_dir)
+ print(f"Opt: {Edot = }")
+
+ def test_Edot_opt_6(self):
+ analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\analysis_temp\ga_airfoil_7"
+ Edot = calculate_CPK_power_consumption(analysis_dir)
+ print(f"{Edot = }")
diff --git a/pymead/tests/aero_tests/mechanical_power_test.py b/pymead/tests/aero_tests/mechanical_power_test.py
index 9d9fb139..41a28ee2 100644
--- a/pymead/tests/aero_tests/mechanical_power_test.py
+++ b/pymead/tests/aero_tests/mechanical_power_test.py
@@ -73,6 +73,10 @@ def test_PK_calculation(self):
plt.show()
pass
+ def test_CPK_calculation_baseline(self):
+ analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_70\analysis\analysis_0"
+ calculate_CPK_mses(analysis_dir)
+
def test_PK_calculation_opt_pai(self):
analysis_dir = r"C:\Users\mlauer2\Documents\pymead\pymead\pymead\tests\pai\root_underwing_opt\opt_runs\2023_05_03_A\ga_opt_61\analysis\analysis_284_w0-100"
# fig, ax = plt.subplots()