diff --git a/src/VSM/core/BodyAerodynamics.py b/src/VSM/core/BodyAerodynamics.py index 4eaf8ac..a4d190e 100644 --- a/src/VSM/core/BodyAerodynamics.py +++ b/src/VSM/core/BodyAerodynamics.py @@ -17,6 +17,7 @@ thin_airfoil_strips, ) from . import jit_cross, jit_norm, jit_dot +from VSM.core import panel_kernels def _rigid_rotation_between(points_0, points_1, dirs_0, dirs_1, rel_tol=1e-9): @@ -1298,17 +1299,17 @@ def compute_trefftz_plane_induced_drag(self, gamma, rho): ).astype(float) d_gamma = np.concatenate(([gamma[0]], gamma[1:] - gamma[:-1], [-gamma[-1]])) q = origins - np.outer(origins @ e_inf, e_inf) # projected onto the plane - drag = 0.0 - for i in range(n): - mid = 0.5 * (q[i] + q[i + 1]) - d = mid - q - r2 = np.sum(d * d, axis=1) - v_t = np.sum( - (d_gamma / (2.0 * np.pi * r2))[:, None] * np.cross(e_inf, d), axis=0 - ) - # trace segment in the sense of the bound vortex (section i+1 -> i) - drag += gamma[i] * np.dot(np.cross(v_t, q[i] - q[i + 1]), e_inf) - return 0.5 * rho * drag + # velocity of all wake vortices k at every trace-segment midpoint i + mid = 0.5 * (q[:-1] + q[1:]) # (n, 3) + d = mid[:, None, :] - q[None, :, :] # (n, n+1, 3) + r2 = np.sum(d * d, axis=2) + v_t = np.sum( + (d_gamma[None, :] / (2.0 * np.pi * r2))[:, :, None] * np.cross(e_inf, d), + axis=1, + ) + # trace segment in the sense of the bound vortex (section i+1 -> i) + drag = np.sum(gamma[:n] * (np.cross(v_t, q[:-1] - q[1:]) @ e_inf)) + return 0.5 * rho * float(drag) def compute_attached_trailed_vortex_forces( self, gamma, rho, core_radius_fraction, va_array @@ -1331,11 +1332,7 @@ def compute_attached_trailed_vortex_forces( (forces, points): ``(n+1, 3)`` force vectors and the points they act at. """ panels = self.panels - n = len(panels) gamma = np.asarray(gamma, dtype=float).ravel() - gamma_net = np.zeros(n + 1) - gamma_net[:n] += gamma # leg bound_1 -> TE_1 of panel j sits on boundary j - gamma_net[1:] -= gamma # leg TE_2 -> bound_2 of panel j sits on boundary j+1 bound_point_1 = np.ascontiguousarray( [p.bound_point_1 for p in panels], dtype=float @@ -1355,13 +1352,7 @@ def compute_attached_trailed_vortex_forces( t = (cp - ac) / (1.0 - ac) points = np.ascontiguousarray(bound + t * legs) - # inflow at a boundary: mean of the neighbouring panels' inflow va_array = np.asarray(va_array, dtype=float) - va_boundary = np.zeros((n + 1, 3)) - va_boundary[:n] += va_array - va_boundary[1:] += va_array - va_boundary[1:n] *= 0.5 - va_speeds = np.linalg.norm(va_array, axis=1) wake_units, wake_speeds = self._wake_directions( va_speeds, va_array / va_speeds[:, None] @@ -1377,8 +1368,9 @@ def compute_attached_trailed_vortex_forces( wake_speeds, float(core_radius_fraction), ) - v_rel = va_boundary + v_ind - forces = rho * gamma_net[:, None] * np.cross(v_rel, legs) + # net leg circulation and boundary inflow (mean of the neighbouring + # panels') are formed in the kernel + forces = panel_kernels.attached_trailed_forces(gamma, va_array, v_ind, legs, rho) return forces, points def compute_results( @@ -1416,10 +1408,6 @@ def compute_results( cl_array[icp] = panel_i.compute_cl(alpha_array[icp]) cd_array[icp], cm_array[icp] = panel_i.compute_cd_cm(alpha_array[icp]) panel_width_array[icp] = panel_i.width - lift = (cl_array * 0.5 * rho * Umag_array**2 * chord_array)[:, np.newaxis] - drag = (cd_array * 0.5 * rho * Umag_array**2 * chord_array)[:, np.newaxis] - moment = (cm_array * 0.5 * rho * Umag_array**2 * chord_array**2)[:, np.newaxis] - # Quarter-chord force directions (TAT3) apply to the VSM model, whose # circulation loop evaluates at 3/4 chord. The LLT model already # evaluates on the quarter chord without any bound-vortex influence @@ -1440,6 +1428,28 @@ def compute_results( alpha_corrected = alpha_array[:, np.newaxis] alpha_uncorrected = alpha_array[:, np.newaxis] + # Section force law (VSM.core.panel_kernels): magnitudes from the polar + # at the 3/4-chord angle and section speed, directions from + # alpha_corrected (the quarter-chord angle, or the 3/4-chord one). + lift_2D_array, drag_2D_array, moment_2D_array, dir_flow_array = ( + panel_kernels.section_force_vectors( + alpha_array, + Umag_array, + cl_array, + cd_array, + cm_array, + chord_array, + x_airf_array, + y_airf_array, + z_airf_array, + rho, + alpha_dir=np.ravel(alpha_corrected), + ) + ) + moment_local_3D_array = panel_kernels.panel_moment_local( + moment_2D_array, panel_width_array, z_airf_array + ) + # Initializing variables cl_prescribed_va_list = [] cd_prescribed_va_list = [] @@ -1496,45 +1506,19 @@ def compute_results( ### Defining panel_variables # Defining directions of airfoil that defines current panel_i z_airf_span = panel_i.z_airf # along the span - y_airf_chord = panel_i.y_airf # along the chord - x_airf_normal_to_chord = panel_i.x_airf # normal to the chord - # TODO: implement these - alpha_corrected_i = alpha_corrected[i] panel_chord = panel_i.chord panel_width = panel_i.width panel_area = panel_chord * panel_width area_all_panels += panel_area - ### Calculate the direction of the induced apparent wind speed to the airfoil orientation - # this is done using the CORRECTED CALCULATED (comes from gamma distribution) angle of attack - # For VSM the correction is applied, and it is the angle of attack, from calculating induced velocities at the 1/4c aerodynamic center location - # For LTT the correction is NOT applied, and it is the angle of attack, from calculating induced velocities at the 3/4c control point - induced_va_airfoil = ( - np.cos(alpha_corrected_i) * y_airf_chord - + np.sin(alpha_corrected_i) * x_airf_normal_to_chord - ) - dir_induced_va_airfoil = induced_va_airfoil / jit_norm(induced_va_airfoil) - - ### Calculate the direction of the lift and drag vectors - # lift is perpendical/normal to induced apparent wind speed - dir_lift_induced_va = jit_cross(dir_induced_va_airfoil, z_airf_span) - dir_lift_induced_va = dir_lift_induced_va / jit_norm(dir_lift_induced_va) - # drag is parallel/tangential to induced apparent wind speed - # (TAT3: the local inner flow direction, which lies in the panel's - # span-perpendicular plane; the global span axis is NOT used here) - dir_drag_induced_va = dir_induced_va_airfoil - - ### Calculating the MAGNITUDE of the lift and drag - # The VSM and LTT methods do NOT differ here, both use the uncorrected angle of attack - # i.e. evaluate the magnitude at the (3/4c) control point - # 2D AIRFOIL aerodynamic forces, so multiplied by chord - lift_induced_va_mag = lift[i] - drag_induced_va_mag = drag[i] - - # panel force VECTOR NORMAL to CALCULATED induced velocity - lift_induced_va = lift_induced_va_mag * dir_lift_induced_va - # panel force VECTOR TANGENTIAL to CALCULATED induced velocity - drag_induced_va = drag_induced_va_mag * dir_drag_induced_va + ### Section force vectors (per unit span) from the force law above: + # lift normal and drag along the direction-defining flow + # dir_induced_va_airfoil (TAT3: the local inner flow direction, in + # the panel's span-perpendicular plane; the global span axis is NOT + # used), magnitudes at the 3/4c control point for VSM and LLT alike. + dir_induced_va_airfoil = dir_flow_array[i] + lift_induced_va = lift_2D_array[i] + drag_induced_va = drag_2D_array[i] ftotal_induced_va = lift_induced_va + drag_induced_va ### Converting forces to prescribed wing va @@ -1689,13 +1673,10 @@ def compute_results( # (1) Panel aerodynamic center in global frame: panel_ac_global = panel_i.aerodynamic_center # 3D [x, y, z] - # (2) Convert local (2D) pitching moment to a 3D vector in global coords. - # Use the axis around which the moment is defined, - # which is the z-axis pointing "spanwise" - moment_axis_global = panel_i.z_airf - - # Scale by panel width if your 'moment[i]' is 2D moment-per-unit-span: - M_local_3D = moment[i] * moment_axis_global * panel_i.width + # (2) Local (2D) pitching moment as a 3D vector in global coords: + # along the z-axis pointing "spanwise", times the width + # (panel_kernels.panel_moment_local). + M_local_3D = moment_local_3D_array[i] # (3) Force application point used for r x F: # use aerodynamic center (1/4c) for all models. @@ -1937,7 +1918,7 @@ def compute_results( ( "cm_panel_dist", np.divide( - moment[:, 0], + moment_2D_array, 0.5 * rho * Umag_array**2 * chord_array**2, out=np.zeros_like(chord_array), where=(Umag_array > 0.0) & (chord_array > 0.0), @@ -2144,39 +2125,39 @@ def _curved_bound_correction_cached( } return AIC_corr - def update_effective_angle_of_attack_if_VSM( - self, - gamma, - core_radius_fraction, - x_airf_array, - y_airf_array, - va_array, - va_norm_array, - va_unit_array, + def compute_quarter_chord_AIC_matrices( + self, core_radius_fraction, va_norm_array, va_unit_array ): - """Updates the angle of attack at the aerodynamic center of each panel, - Calculated at the AERODYNAMIC CENTER, which needs an update for VSM - And can just use the old value for the LLT + """Induction matrices of the force-direction flow at the panel + quarter chords (aerodynamic centres): ``v_qc = va + AIC_qc gamma``. + + The direction-defining flow is the undisturbed relative flow at the + quarter chord (TAT3 in Gaunaa, Li & Pirrung, TORQUE 2026). The LLT + AIC holds the freestream-independent trailed system (legs + wake) + evaluated on the quarter chord, which is Prandtl's exact downwash for + a straight line but misses the influence of a curved or swept bound + vortex. Evaluating the other panels' bound segments on the line itself + is logarithmically singular (Li et al. 2020, "LL-P&S"), so the + curvature/sweep/skew influence is added as the chord-regularised + difference between the actual and the locally unrolled vortex + systems; see assemble_curved_bound_correction_AIC. + + The matrices depend on the geometry and the wake (``va_norm_array``, + ``va_unit_array``, as in :meth:`compute_AIC_matrices`), not on the + circulation, so a caller can freeze them numerically, e.g. as + parameters of a CasADi graph (see ``VSM.core.panel_kernels``). Args: - None + core_radius_fraction (float): Vortex core radius fraction. + va_norm_array (np.ndarray): (n,) panel inflow speeds [m/s]. + va_unit_array (np.ndarray): (n, 3) panel inflow directions. Returns: - None + Tuple[np.ndarray, np.ndarray, np.ndarray]: the x, y and z + components, each (n, n): ``AIC_LLT + AIC_curved_correction``. """ - # The direction-defining flow is the undisturbed relative flow at the - # quarter chord (TAT3 in Gaunaa, Li & Pirrung, TORQUE 2026). The LLT - # AIC holds the freestream-independent trailed system (legs + wake) - # evaluated on the quarter chord, which is Prandtl's exact downwash for - # a straight line but misses the influence of a curved or swept bound - # vortex. Evaluating the other panels' bound segments on the line itself - # is logarithmically singular (Li et al. 2020, "LL-P&S"), so the - # curvature/sweep/skew influence is added as the chord-regularised - # difference between the actual and the locally unrolled vortex - # systems; see assemble_curved_bound_correction_AIC. - aerodynamic_model_type = "LLT" AIC_x, AIC_y, AIC_z = self.compute_AIC_matrices( - aerodynamic_model_type, core_radius_fraction, va_norm_array, va_unit_array + "LLT", core_radius_fraction, va_norm_array, va_unit_array ) eval_points = np.ascontiguousarray( [panel.aerodynamic_center for panel in self.panels], dtype=float @@ -2194,7 +2175,9 @@ def update_effective_angle_of_attack_if_VSM( [panel.TE_point_2 for panel in self.panels], dtype=float ) chords = np.ascontiguousarray([panel.chord for panel in self.panels], dtype=float) - chord_dir = np.ascontiguousarray(y_airf_array, dtype=float) + chord_dir = np.ascontiguousarray( + [panel.y_airf for panel in self.panels], dtype=float + ) # The correction depends on the panel geometry only (the wake terms # cancel between the actual and the unrolled systems, and the leg core # radius has no measurable effect), so it is computed once per shape @@ -2209,21 +2192,34 @@ def update_effective_angle_of_attack_if_VSM( chord_dir, chords, np.ascontiguousarray(va_norm_array, dtype=float).ravel(), float(core_radius_fraction), ) - AIC_x = AIC_x + AIC_corr[0] - AIC_y = AIC_y + AIC_corr[1] - AIC_z = AIC_z + AIC_corr[2] - induced_velocity_all = np.array( - [ - np.matmul(AIC_x, gamma), - np.matmul(AIC_y, gamma), - np.matmul(AIC_z, gamma), - ] - ).T - relative_velocity_array = va_array + induced_velocity_all - v_normal_array = np.sum(x_airf_array * relative_velocity_array, axis=1) - v_tangential_array = np.sum(y_airf_array * relative_velocity_array, axis=1) - alpha_array = np.arctan2(v_normal_array, v_tangential_array) + return AIC_x + AIC_corr[0], AIC_y + AIC_corr[1], AIC_z + AIC_corr[2] + + def update_effective_angle_of_attack_if_VSM( + self, + gamma, + core_radius_fraction, + x_airf_array, + y_airf_array, + va_array, + va_norm_array, + va_unit_array, + ): + """Force-direction angle of each panel: the angle of the relative flow + at the quarter chord, ``va + AIC_qc gamma`` with the matrices of + :meth:`compute_quarter_chord_AIC_matrices`. + Returns: + np.ndarray: (n, 1) angles [rad]. + """ + AIC_x, AIC_y, AIC_z = self.compute_quarter_chord_AIC_matrices( + core_radius_fraction, va_norm_array, va_unit_array + ) + relative_velocity_array = panel_kernels.quarter_chord_velocity( + va_array, gamma, AIC_x, AIC_y, AIC_z + ) + alpha_array = panel_kernels.direction_angle( + relative_velocity_array, x_airf_array, y_airf_array + ) return alpha_array[:, np.newaxis] def compute_line_aerodynamic_force( diff --git a/src/VSM/core/Solver.py b/src/VSM/core/Solver.py index 576f3ba..96d3766 100644 --- a/src/VSM/core/Solver.py +++ b/src/VSM/core/Solver.py @@ -3,6 +3,7 @@ import logging from scipy.linalg import solve_banded from . import jit_cross +from VSM.core import panel_kernels class Solver: @@ -201,6 +202,7 @@ def __init__( # rejected and the step quartered (0.5 = residual may at most double). self._newton_reject_ratio = 0.5 self._casadi_newton_cache: dict = {} + self._casadi_newton_loop_cache: dict = {} #: Diagnostic: did the last casadi_newton solve fall back to base? self.last_fallback = False #: Diagnostic: residual evaluations used by the last casadi_newton solve. @@ -1028,56 +1030,12 @@ def _casadi_newton_key(self, alpha_tables, cl_tables, cd_tables, cm_tables) -> t digest.hexdigest(), ) - def _build_casadi_newton_function(self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None): - """Build the CasADi function of the PER-PANEL section physics that - :meth:`gamma_loop_casadi_newton` assembles its residual and exact - Jacobian from. - - The circulation residual whose root the Newton loop finds is - - R(gamma) = (I - diag(mu(alpha)) L) gamma - G_raw(gamma), - - exactly the fixed point the ``base`` and ``anderson`` loops iterate - towards (``gamma = (I - diag(mu) L)^-1 G_raw(gamma)``, see - :meth:`_regularize_gamma_target`), including the Li/Gaunaa viscosity - ``mu_i = max(0, -k S Cl'_i / dz_i^2)`` gated on any panel being past - its stall onset (``mu = 0`` without artificial viscosity). Everything - that is not the dense induction ``v_ind = AIC gamma`` is local to a - panel: ``G_raw_i``, ``mu_i`` and ``alpha_i`` depend on ``gamma`` only - through the panel's own relative velocity ``v_rel_i``. So the symbolic - function takes ``v_rel`` (and ``L gamma``) as INPUTS and returns - - h_i = G_raw_i + mu_i (L gamma)_i, c_i = dh_i / dv_rel_i (3 values) - - plus ``mu``, ``alpha`` and ``|v_rel x z|``; the Newton loop then forms - ``R = gamma - h`` and, by the chain rule through ``v_rel = va + AIC - gamma`` (with ``dh_i/d(L gamma)_i = mu_i``), - - J = I - diag(mu) L - sum_k diag(c_k) AIC_k - - in NumPy. This keeps the CasADi graph O(n) -- independent of the - n x n induction matrices -- so it costs ~0.1 ms per evaluation and - builds in tens of milliseconds, while the dense algebra runs in BLAS. - - The 2-D polars are ``ca.interpolant`` tables. Queries are clamped to - the table range first, so the lookup reproduces ``np.interp`` (constant - beyond both ends) and the lift slope is the same central difference of - the same piecewise-linear table that :meth:`_lift_slope_from_ctx` - evaluates. When every panel shares one alpha grid (the normal outcome - of batch polar generation) the lookup is a single 2-D interpolant over - (alpha, panel index) evaluated exactly at the integer panel nodes; - otherwise one 1-D interpolant per panel. ``polar_interpolation = - "bspline"`` swaps the piecewise-linear table for a cubic spline - through the same nodes (degree 1 along the panel-index axis, so panel - rows never mix). - - Outputs, in order: ``h``, ``dh/dv_rel`` (n x 3, from the corner- - averaged polar), ``mu``, ``alpha``, ``umag``, ``gate``, then for - callers that build FORCES from the same lookups: ``h_avg`` (h with the - corner-averaged Cl -- differentiate this for a kink-tolerant - Jacobian), ``cl``, ``cd``, ``cm`` at ``alpha`` and their corner- - averaged counterparts ``cl_avg``, ``cd_avg``, ``cm_avg``. With - "bspline" polars the averaged outputs equal the exact ones. + def _casadi_polar_lookups(self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None): + """Clamped CasADi lookups ``alpha (n) -> cl, cd, cm (n)`` of the + per-panel polars, shared by :meth:`_build_casadi_newton_function` and + :meth:`_build_casadi_panel_force_function` (see the former for the + interpolation). Returns ``(cl_of, cd_of, cm_of, smooth)``, ``smooth`` + True for "bspline" polars. """ import casadi as ca @@ -1089,16 +1047,6 @@ def _build_casadi_newton_function(self, alpha_tables, cl_tables, cd_tables=None, f"(panel {i} violates this)." ) - v_rel = ca.MX.sym("v_rel", n, 3) - x_airf = ca.MX.sym("x_airf", n, 3) - y_airf = ca.MX.sym("y_airf", n, 3) - z_airf = ca.MX.sym("z_airf", n, 3) - chord = ca.MX.sym("chord", n) - width = ca.MX.sym("width", n) - stall_angles = ca.MX.sym("stall_angles", n) - stall_angles_neg = ca.MX.sym("stall_angles_neg", n) - l_gamma = ca.MX.sym("l_gamma", n) # (L gamma), a parameter here - shared_grid = all( table.shape == alpha_tables[0].shape and np.array_equal(table, alpha_tables[0]) @@ -1154,6 +1102,75 @@ def lookup(alpha_vec): cl_of = make_lookup("cl", cl_tables) cd_of = make_lookup("cd", cd_tables) cm_of = make_lookup("cm", cm_tables) + return cl_of, cd_of, cm_of, smooth + + def _build_casadi_newton_function(self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None): + """Build the CasADi function of the PER-PANEL section physics that + :meth:`gamma_loop_casadi_newton` assembles its residual and exact + Jacobian from. + + The circulation residual whose root the Newton loop finds is + + R(gamma) = (I - diag(mu(alpha)) L) gamma - G_raw(gamma), + + exactly the fixed point the ``base`` and ``anderson`` loops iterate + towards (``gamma = (I - diag(mu) L)^-1 G_raw(gamma)``, see + :meth:`_regularize_gamma_target`), including the Li/Gaunaa viscosity + ``mu_i = max(0, -k S Cl'_i / dz_i^2)`` gated on any panel being past + its stall onset (``mu = 0`` without artificial viscosity). Everything + that is not the dense induction ``v_ind = AIC gamma`` is local to a + panel: ``G_raw_i``, ``mu_i`` and ``alpha_i`` depend on ``gamma`` only + through the panel's own relative velocity ``v_rel_i``. So the symbolic + function takes ``v_rel`` (and ``L gamma``) as INPUTS and returns + + h_i = G_raw_i + mu_i (L gamma)_i, c_i = dh_i / dv_rel_i (3 values) + + plus ``mu``, ``alpha`` and ``|v_rel x z|``; the Newton loop then forms + ``R = gamma - h`` and, by the chain rule through ``v_rel = va + AIC + gamma`` (with ``dh_i/d(L gamma)_i = mu_i``), + + J = I - diag(mu) L - sum_k diag(c_k) AIC_k + + in NumPy. This keeps the CasADi graph O(n) -- independent of the + n x n induction matrices -- so it costs ~0.1 ms per evaluation and + builds in tens of milliseconds, while the dense algebra runs in BLAS. + + The 2-D polars are ``ca.interpolant`` tables. Queries are clamped to + the table range first, so the lookup reproduces ``np.interp`` (constant + beyond both ends) and the lift slope is the same central difference of + the same piecewise-linear table that :meth:`_lift_slope_from_ctx` + evaluates. When every panel shares one alpha grid (the normal outcome + of batch polar generation) the lookup is a single 2-D interpolant over + (alpha, panel index) evaluated exactly at the integer panel nodes; + otherwise one 1-D interpolant per panel. ``polar_interpolation = + "bspline"`` swaps the piecewise-linear table for a cubic spline + through the same nodes (degree 1 along the panel-index axis, so panel + rows never mix). + + Outputs, in order: ``h``, ``dh/dv_rel`` (n x 3, from the corner- + averaged polar), ``mu``, ``alpha``, ``umag``, ``gate``, then for + callers that build FORCES from the same lookups: ``h_avg`` (h with the + corner-averaged Cl -- differentiate this for a kink-tolerant + Jacobian), ``cl``, ``cd``, ``cm`` at ``alpha`` and their corner- + averaged counterparts ``cl_avg``, ``cd_avg``, ``cm_avg``. With + "bspline" polars the averaged outputs equal the exact ones. + """ + import casadi as ca + + n = self.n_panels + cl_of, cd_of, cm_of, smooth = self._casadi_polar_lookups( + alpha_tables, cl_tables, cd_tables, cm_tables + ) + + v_rel = ca.MX.sym("v_rel", n, 3) + x_airf = ca.MX.sym("x_airf", n, 3) + y_airf = ca.MX.sym("y_airf", n, 3) + z_airf = ca.MX.sym("z_airf", n, 3) + chord = ca.MX.sym("chord", n) + width = ca.MX.sym("width", n) + stall_angles = ca.MX.sym("stall_angles", n) + stall_angles_neg = ca.MX.sym("stall_angles_neg", n) + l_gamma = ca.MX.sym("l_gamma", n) # (L gamma), a parameter here def cross_rows(a, b): return ca.horzcat( @@ -1264,6 +1281,110 @@ def _casadi_newton_function(self): self._casadi_newton_cache[key] = fn return fn + def _build_casadi_panel_force_function( + self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None + ): + """CasADi function of the per-panel section forces, the force law of + ``compute_results`` (``VSM.core.panel_kernels.panel_forces``) on the + polar lookups of the circulation solve. See + :meth:`casadi_panel_force_function`.""" + import casadi as ca + + n = self.n_panels + cl_of, cd_of, cm_of, _ = self._casadi_polar_lookups( + alpha_tables, cl_tables, cd_tables, cm_tables + ) + v_rel = ca.MX.sym("v_rel", n, 3) + v_qc = ca.MX.sym("v_qc", n, 3) + x_airf = ca.MX.sym("x_airf", n, 3) + y_airf = ca.MX.sym("y_airf", n, 3) + z_airf = ca.MX.sym("z_airf", n, 3) + chord = ca.MX.sym("chord", n) + width = ca.MX.sym("width", n) + rho = ca.MX.sym("rho") + + alpha, umag = panel_kernels.relative_flow(v_rel, x_airf, y_airf, z_airf, xp=ca) + alpha_dir = panel_kernels.direction_angle(v_qc, x_airf, y_airf, xp=ca) + cl, cd, cm = cl_of(alpha), cd_of(alpha), cm_of(alpha) + force, moment_local = panel_kernels.panel_forces( + alpha, umag, cl, cd, cm, chord, width, x_airf, y_airf, z_airf, rho, + alpha_dir=alpha_dir, xp=ca, + ) + return ca.Function( + "vsm_panel_forces", + [v_rel, v_qc, x_airf, y_airf, z_airf, chord, width, rho], + [force, moment_local, alpha, alpha_dir, umag, cl, cd, cm], + ["v_rel", "v_qc", "x_airf", "y_airf", "z_airf", "chord", "width", "rho"], + ["force", "moment_local", "alpha", "alpha_dir", "umag", "cl", "cd", "cm"], + ) + + def casadi_panel_force_function(self): + """CasADi function of the per-panel section forces of the current + panels (the ones of the last :meth:`solve`), cached per polar set like + the ``casadi_newton`` residual function. + + Inputs (numeric or symbolic), ``n = n_panels``: + + * ``v_rel`` (n, 3): relative velocity at the 3/4-chord control points, + ``va + AIC gamma`` -- sets ``alpha``, ``umag`` and so the polar + magnitudes ``cl, cd, cm``; + * ``v_qc`` (n, 3): the force-direction flow. For the quarter-chord law + (``is_aoa_corrected=True``) ``va + AIC_qc gamma`` with + ``BodyAerodynamics.compute_quarter_chord_AIC_matrices`` (or + ``panel_kernels.quarter_chord_velocity``); for ``False`` pass + ``v_rel`` again; + * ``x_airf, y_airf, z_airf`` (n, 3), ``chord, width`` (n), ``rho``. + + Outputs: ``force`` (n, 3) [N], ``moment_local`` (n, 3) [Nm] (the + quarter-chord pitching moment along ``z_airf``), ``alpha``, + ``alpha_dir`` [rad], ``umag`` [m/s], ``cl, cd, cm`` (n). + + Same as ``compute_results`` per panel WITHOUT the attached-trailed + vortex force (``panel_kernels.attached_trailed_forces``) and the + viscous drag correction. Taking ``v_qc`` rather than the matrices keeps + the graph O(n); the caller forms the dense product. + """ + import casadi # noqa: F401 (clear error when the optional dep is missing) + + alpha_tables, cl_tables, cd_tables, cm_tables = self._polar_tables() + key = ("panel_forces",) + self._casadi_newton_key( + alpha_tables, cl_tables, cd_tables, cm_tables + ) + fn = self._casadi_newton_cache.get(key) + if fn is None: + fn = self._build_casadi_panel_force_function( + alpha_tables, cl_tables, cd_tables, cm_tables + ) + if len(self._casadi_newton_cache) >= 32: + self._casadi_newton_cache.pop(next(iter(self._casadi_newton_cache))) + self._casadi_newton_cache[key] = fn + return fn + + def _casadi_newton_loop_function(self): + """The residual function restricted to the six outputs the Newton loop + reads (``h, dh_dvrel, mu, alpha, umag, gate``): the same graph, minus + the force-side lookups it would otherwise evaluate on every call. + Cached per polar set alongside the full function.""" + import casadi as ca + + alpha_tables, cl_tables, cd_tables, cm_tables = self._polar_tables() + key = self._casadi_newton_key(alpha_tables, cl_tables, cd_tables, cm_tables) + loop_fn = self._casadi_newton_loop_cache.get(key) + if loop_fn is None: + fn = self._casadi_newton_function() + inputs = fn.mx_in() + loop_fn = ca.Function( + "vsm_panel_residual_loop", + inputs, + fn.call(inputs)[:6], + fn.name_in(), + fn.name_out()[:6], + ) + if len(self._casadi_newton_loop_cache) >= 32: + self._casadi_newton_loop_cache.pop(next(iter(self._casadi_newton_loop_cache))) + self._casadi_newton_loop_cache[key] = loop_fn + return loop_fn + def gamma_loop_casadi_newton(self, gamma_initial: np.ndarray) -> tuple: """Newton solve of the circulation residual with an exact CasADi Jacobian, globalised by pseudo-transient continuation @@ -1312,7 +1433,9 @@ def gamma_loop_casadi_newton(self, gamma_initial: np.ndarray) -> tuple: fallback to decide. """ try: - fn = self._casadi_newton_function() + import casadi as ca + + fn = self._casadi_newton_loop_function() except ImportError as exc: # pragma: no cover - depends on optional dep raise ImportError( "gamma_loop_type='casadi_newton' requires the optional 'casadi' " @@ -1329,14 +1452,19 @@ def gamma_loop_casadi_newton(self, gamma_initial: np.ndarray) -> tuple: stall_angles = np.full(n, np.inf) stall_angles_neg = np.full(n, -np.inf) laplacian = np.zeros((n, n)) - params = ( - self.x_airf_array, - self.y_airf_array, - self.z_airf_array, - self.chord_array, - self.width_array, - stall_angles, - stall_angles_neg, + # Constant within the solve: converted to CasADi matrices once, which + # halves the cost of each residual evaluation. + params = tuple( + ca.DM(np.asarray(value, dtype=float)) + for value in ( + self.x_airf_array, + self.y_airf_array, + self.z_airf_array, + self.chord_array, + self.width_array, + stall_angles, + stall_angles_neg, + ) ) aic = (self.AIC_x, self.AIC_y, self.AIC_z) @@ -1345,7 +1473,7 @@ def evaluate(x): # v_rel = va + AIC gamma (see _build_casadi_newton_function). v_rel = self.va_array + np.column_stack([a @ x for a in aic]) l_gamma = laplacian @ x - h, dh_dvrel, mu, alpha, umag, gate = fn(v_rel, *params, l_gamma)[:6] + h, dh_dvrel, mu, alpha, umag, gate = fn(v_rel, *params, l_gamma) h = np.asarray(h, dtype=float).ravel() dh_dvrel = np.asarray(dh_dvrel, dtype=float) mu = np.asarray(mu, dtype=float).ravel() diff --git a/src/VSM/core/WingGeometry.py b/src/VSM/core/WingGeometry.py index ecc56c4..81fb1fd 100644 --- a/src/VSM/core/WingGeometry.py +++ b/src/VSM/core/WingGeometry.py @@ -627,22 +627,18 @@ def compute_projected_area(self, z_plane_vector=np.array([0.0, 0.0, 1.0])): """ sections = self.sections n = np.asarray(z_plane_vector, dtype=float) - n /= np.linalg.norm(n) - - def proj(P): - # orthogonal projection onto plane with normal n - return P - np.dot(P, n) * n - - S = 0.0 - for i in range(len(sections) - 1): - A = proj(sections[i].LE_point) - B = proj(sections[i].TE_point) - C = proj(sections[i + 1].TE_point) - D = proj(sections[i + 1].LE_point) - # two triangles: (A,B,C) and (A,C,D) - S += 0.5 * np.linalg.norm(np.cross(B - A, C - A)) - S += 0.5 * np.linalg.norm(np.cross(C - A, D - A)) - return S + n = n / np.linalg.norm(n) # not in place: never modify the caller's array + + LE = np.array([section.LE_point for section in sections], dtype=float) + TE = np.array([section.TE_point for section in sections], dtype=float) + # orthogonal projection onto the plane with normal n + LE = LE - np.outer(LE @ n, n) + TE = TE - np.outer(TE @ n, n) + A, B, C, D = LE[:-1], TE[:-1], TE[1:], LE[1:] + # two triangles per strip: (A,B,C) and (A,C,D) + area_ABC = 0.5 * np.linalg.norm(np.cross(B - A, C - A), axis=1) + area_ACD = 0.5 * np.linalg.norm(np.cross(C - A, D - A), axis=1) + return float(np.sum(np.column_stack([area_ABC, area_ACD]))) def compute_flat_area(self) -> float: """Calculate the total flat (undeformed) wing area in 3D space. diff --git a/src/VSM/core/panel_kernels.py b/src/VSM/core/panel_kernels.py new file mode 100644 index 0000000..f074f61 --- /dev/null +++ b/src/VSM/core/panel_kernels.py @@ -0,0 +1,236 @@ +"""Section force law of the VSM as ``xp`` kernels. + +Pure functions over a math namespace ``xp``, either ``numpy`` (the default) +or ``casadi``, so the same formula serves ``BodyAerodynamics.compute_results`` +and a CasADi graph (``Solver.casadi_panel_force_function``, or a caller's own +graph such as AWETrim's CasADi trim). The formula lives here only; callers do +not restate it. + +Conventions: + +* Per-panel vectors are rows of an ``(n, 3)`` array; per-panel scalars are + ``(n,)`` arrays in NumPy and ``(n, 1)`` columns in CasADi. +* ``x_airf`` is the chord normal, ``y_airf`` the chord (leading to trailing + edge), ``z_airf`` the span direction of each panel (``Panel`` axes). +* Angles are in radians. ``alpha = atan2(v . x_airf, v . y_airf)``. +* ``umag = |v_rel x z_airf|`` is the speed the 2-D section sees; the dynamic + pressure of the section force is ``0.5 rho umag^2``, not the freestream one. + +The force law (Gaunaa, Li & Pirrung, TORQUE 2026): the MAGNITUDES of lift, +drag and moment come from the polar at the 3/4-chord control-point angle +``alpha``; the DIRECTIONS come from the direction angle ``alpha_dir``. With +``is_aoa_corrected=True`` (TAT3, VSM model) ``alpha_dir`` is the angle of the +flow at the quarter chord, ``va + AIC_qc gamma`` with the matrices of +``BodyAerodynamics.compute_quarter_chord_AIC_matrices``; with ``False`` it is +``alpha`` itself. +""" + +from __future__ import annotations + +from typing import Any + +import numpy as np + + +# --------------------------------------------------------------------- xp rows +def _is_casadi(xp: Any) -> bool: + return hasattr(xp, "horzcat") + + +def hstack_columns(xp: Any, *cols): + """``(n, k)`` from ``k`` per-row columns.""" + if _is_casadi(xp): + return xp.horzcat(*cols) + return np.column_stack(cols) + + +def cross_rows(a, b, xp=np): + """Row-wise cross product of two ``(n, 3)`` arrays.""" + return hstack_columns( + xp, + a[:, 1] * b[:, 2] - a[:, 2] * b[:, 1], + a[:, 2] * b[:, 0] - a[:, 0] * b[:, 2], + a[:, 0] * b[:, 1] - a[:, 1] * b[:, 0], + ) + + +def dot_rows(a, b, xp=np): + """Row-wise dot product of two ``(n, 3)`` arrays -> per-row scalar.""" + if _is_casadi(xp): + return xp.sum2(a * b) + return np.sum(a * b, axis=1) + + +def norm_rows(a, xp=np): + return xp.sqrt(dot_rows(a, a, xp)) + + +def scale_rows(s, a, xp=np): + """Per-row scalar ``s`` times the ``(n, 3)`` rows of ``a``.""" + if _is_casadi(xp): + return xp.repmat(s, 1, a.shape[1]) * a + return np.asarray(s)[:, None] * a + + +def unit_rows(a, xp=np): + """Rows of ``a`` divided by their norms.""" + norms = norm_rows(a, xp) + if _is_casadi(xp): + return a / xp.repmat(norms, 1, a.shape[1]) + return a / norms[:, None] + + +def _atan2(y, x, xp): + return xp.atan2(y, x) if _is_casadi(xp) else np.arctan2(y, x) + + +def _matvec(matrix, vector, xp): + if _is_casadi(xp): + return xp.mtimes(matrix, vector) + return matrix @ vector + + +# ------------------------------------------------------------------ section flow +def relative_flow(v_rel, x_airf, y_airf, z_airf, xp=np): + """Angle of attack and section speed of the relative velocity ``v_rel`` + (``Solver.compute_aerodynamic_quantities``): + ``alpha = atan2(v_rel . x_airf, v_rel . y_airf)``, ``umag = |v_rel x z_airf|``.""" + alpha = direction_angle(v_rel, x_airf, y_airf, xp) + umag = norm_rows(cross_rows(v_rel, z_airf, xp), xp) + return alpha, umag + + +def quarter_chord_velocity(va, gamma, AIC_x, AIC_y, AIC_z, xp=np): + """Relative flow at the panel quarter chords, ``va + AIC_qc gamma``. + + ``AIC_*`` are the ``(n, n)`` quarter-chord induction matrices of + ``BodyAerodynamics.compute_quarter_chord_AIC_matrices`` (trailed system on + the quarter chord plus the curved bound-vortex correction), ``va`` the + ``(n, 3)`` panel inflow and ``gamma`` the ``(n,)`` circulation. + """ + induced = hstack_columns( + xp, + _matvec(AIC_x, gamma, xp), + _matvec(AIC_y, gamma, xp), + _matvec(AIC_z, gamma, xp), + ) + return va + induced + + +def direction_angle(v, x_airf, y_airf, xp=np): + """Angle of the flow ``v`` (``(n, 3)``) in each panel's section plane, + ``atan2(v . x_airf, v . y_airf)``. Applied to the quarter-chord flow it is + the TAT3 direction angle of the force law.""" + return _atan2(dot_rows(x_airf, v, xp), dot_rows(y_airf, v, xp), xp) + + +# ------------------------------------------------------------------- force law +def section_force_directions(alpha_dir, x_airf, y_airf, z_airf, xp=np): + """Unit flow and lift directions of each section at direction angle + ``alpha_dir``: ``dir_flow = unit(cos(a) y_airf + sin(a) x_airf)`` (in the + plane normal to ``z_airf``; drag acts along it, TAT3's local inner flow, + not the global span axis) and ``dir_lift = unit(dir_flow x z_airf)``.""" + dir_flow = scale_rows(xp.cos(alpha_dir), y_airf, xp) + scale_rows( + xp.sin(alpha_dir), x_airf, xp + ) + dir_flow = unit_rows(dir_flow, xp) + dir_lift = unit_rows(cross_rows(dir_flow, z_airf, xp), xp) + return dir_flow, dir_lift + + +def section_force_magnitudes(umag, cl, cd, cm, chord, rho): + """2-D lift, drag [N/m] and pitching moment [Nm/m] about the quarter + chord, ``c_x 0.5 rho umag^2 chord`` (``chord^2`` for the moment). The + coefficients are the polar values at the 3/4-chord angle.""" + lift = cl * 0.5 * rho * umag**2 * chord + drag = cd * 0.5 * rho * umag**2 * chord + moment = cm * 0.5 * rho * umag**2 * chord**2 + return lift, drag, moment + + +def section_force_vectors( + alpha, umag, cl, cd, cm, chord, x_airf, y_airf, z_airf, rho, alpha_dir=None, xp=np +): + """2-D (per unit span) lift and drag vectors and the flow direction. + + ``alpha_dir`` None uses ``alpha`` (the ``is_aoa_corrected=False`` law). + Returns ``(lift_2d (n, 3), drag_2d (n, 3), moment (n,), dir_flow (n, 3))``. + """ + if alpha_dir is None: + alpha_dir = alpha + lift, drag, moment = section_force_magnitudes(umag, cl, cd, cm, chord, rho) + dir_flow, dir_lift = section_force_directions(alpha_dir, x_airf, y_airf, z_airf, xp) + return scale_rows(lift, dir_lift, xp), scale_rows(drag, dir_flow, xp), moment, dir_flow + + +def panel_forces( + alpha, + umag, + cl, + cd, + cm, + chord, + width, + x_airf, + y_airf, + z_airf, + rho, + alpha_dir=None, + xp=np, +): + """Per-panel 3-D force [N] and local pitching moment [Nm]. + + ``force = width (lift_2d + drag_2d)`` (see :func:`section_force_vectors`); + ``moment_local = width (cm q chord^2) z_airf``, the quarter-chord moment as + a vector along the span axis. ``alpha_dir`` None gives the 3/4-chord + directions (``is_aoa_corrected=False``); pass :func:`direction_angle` of + :func:`quarter_chord_velocity` for the quarter-chord law (``True``). + + Not included: the attached-trailed vortex force + (:func:`attached_trailed_forces`) and the viscous drag correction + (``BodyAerodynamics.viscous_drag_correction``, NumPy only). + + Returns ``(force (n, 3), moment_local (n, 3))``. + """ + lift_2d, drag_2d, moment, _ = section_force_vectors( + alpha, umag, cl, cd, cm, chord, x_airf, y_airf, z_airf, rho, alpha_dir, xp + ) + force = scale_rows(width, lift_2d + drag_2d, xp) + return force, panel_moment_local(moment, width, z_airf, xp) + + +def panel_moment_local(moment, width, z_airf, xp=np): + """3-D local pitching moment [Nm] of each panel: the 2-D quarter-chord + moment [Nm/m] (:func:`section_force_magnitudes`) as a vector along the + span axis ``z_airf``, times the panel width. ``(n, 3)``.""" + return scale_rows(width, scale_rows(moment, z_airf, xp), xp) + + +# ---------------------------------------------------- attached trailed vortices +def attached_trailed_forces(gamma, va_panels, v_induced, legs, rho, xp=np): + """Kutta-Joukowski force on the attached trailed (AT) legs (Gaunaa, Li & + Pirrung 2026, Sec. 3), ``rho (gamma_j - gamma_{j-1}) (v_rel_j x l_j)`` on + boundary ``j = 0..n`` with ``gamma_{-1} = gamma_n = 0``. + + ``v_rel_j`` is the mean of the neighbouring panels' inflow ``va_panels`` + (the tip boundaries take their one panel's) plus ``v_induced``, the full + vortex system's induction at the AT force points, ``(n+1, 3)``. ``legs`` + are the ``(n+1, 3)`` leg vectors, bound vortex to trailing edge. Returns + the ``(n+1, 3)`` boundary forces. + """ + n = va_panels.shape[0] + to_boundary = np.zeros((n + 1, n)) # gamma_net = to_boundary @ gamma + to_boundary[np.arange(n), np.arange(n)] += 1.0 + to_boundary[np.arange(1, n + 1), np.arange(n)] -= 1.0 + # boundary inflow = mean of the neighbouring panels (halving is exact) + mean = np.abs(to_boundary) + mean[1:n] *= 0.5 + if _is_casadi(xp): + gamma_net = xp.mtimes(xp.DM(to_boundary), gamma) + va_boundary = xp.mtimes(xp.DM(mean), va_panels) + if isinstance(legs, np.ndarray): + legs = xp.DM(legs) + else: + gamma_net = to_boundary @ np.asarray(gamma, dtype=float).ravel() + va_boundary = mean @ va_panels + return scale_rows(rho * gamma_net, cross_rows(va_boundary + v_induced, legs, xp), xp) diff --git a/src/VSM/core/utils.py b/src/VSM/core/utils.py index d4c4643..10c3630 100644 --- a/src/VSM/core/utils.py +++ b/src/VSM/core/utils.py @@ -30,156 +30,210 @@ def jit_dot(a, b): # round-off; Filament.py remains the readable reference implementation. # The only intentional difference: the bound-vortex in-core branch cannot # emit the logging.info diagnostic from compiled code. +# +# Vectors are carried as 3-tuples of floats, not length-3 arrays: numba +# keeps tuples on the stack, whereas every small-array temporary (np.cross, +# r1 - r0, np.zeros(3)) is a heap allocation, and those allocations were +# ~95% of the kernel cost (~1.1 us -> ~25 ns per filament evaluation). +# error_model="numpy" keeps NumPy's inf/nan semantics for a zero +# denominator (numba's default would raise ZeroDivisionError on scalars), +# so degenerate geometry behaves as it did with the array kernels. # ===================================================================== _ALPHA0_OSEEN = 1.25643 # Oseen parameter (Filament._alpha0) _NU_AIR = 1.48e-5 # kinematic viscosity of air [m^2/s] (Filament._nu) +_ZERO3 = (0.0, 0.0, 0.0) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _row(a, i): + """Row ``i`` of an (n, 3) array as a 3-tuple.""" + return (a[i, 0], a[i, 1], a[i, 2]) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _add(a, b): + return (a[0] + b[0], a[1] + b[1], a[2] + b[2]) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _sub(a, b): + return (a[0] - b[0], a[1] - b[1], a[2] - b[2]) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _scale(s, a): + return (s * a[0], s * a[1], s * a[2]) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _div(a, s): + return (a[0] / s, a[1] / s, a[2] / s) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _dot(a, b): + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2] -@jit(nopython=True, cache=True) + +@jit(nopython=True, cache=True, error_model="numpy") +def _cross(a, b): + return ( + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0], + ) + + +@jit(nopython=True, cache=True, error_model="numpy") +def _norm(a): + return np.sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) + + +@jit(nopython=True, cache=True, error_model="numpy") def _vel_bound_vortex(XV1, XV2, XVP, gamma, core_radius_fraction): """Bound-filament induced velocity (Vatistas core), as in Filament.velocity_3D_bound_vortex.""" - r0 = XV2 - XV1 - r1 = XVP - XV1 - r2 = XVP - XV2 + r0 = _sub(XV2, XV1) + r1 = _sub(XVP, XV1) + r2 = _sub(XVP, XV2) - r1Xr0 = np.cross(r1, r0) - epsilon = core_radius_fraction * np.linalg.norm(r0) - dist = np.linalg.norm(r1Xr0) / np.linalg.norm(r0) + r1Xr0 = _cross(r1, r0) + epsilon = core_radius_fraction * _norm(r0) + dist = _norm(r1Xr0) / _norm(r0) if dist > epsilon: - r1Xr2 = np.cross(r1, r2) - return ( + r1Xr2 = _cross(r1, r2) + return _scale( gamma / (4 * np.pi) - * r1Xr2 - / (np.linalg.norm(r1Xr2) ** 2) - * np.dot(r0, r1 / np.linalg.norm(r1) - r2 / np.linalg.norm(r2)) + / (_norm(r1Xr2) ** 2) + * _dot(r0, _sub(_div(r1, _norm(r1)), _div(r2, _norm(r2)))), + r1Xr2, ) elif dist < 1e-12 * epsilon: - return np.zeros(3) + return _ZERO3 else: - r1_radial = r1 - np.dot(r1, r0) * r0 / (np.linalg.norm(r0) ** 2) - r2_radial = r2 - np.dot(r2, r0) * r0 / (np.linalg.norm(r0) ** 2) - r1_proj = np.dot(r1, r0) * r0 / ( - np.linalg.norm(r0) ** 2 - ) + epsilon * r1_radial / np.linalg.norm(r1_radial) - r2_proj = np.dot(r2, r0) * r0 / ( - np.linalg.norm(r0) ** 2 - ) + epsilon * r2_radial / np.linalg.norm(r2_radial) - r1Xr2_proj = np.cross(r1_proj, r2_proj) - vel_ind_proj = ( + r0_sq = _norm(r0) ** 2 + r1_axial = _div(_scale(_dot(r1, r0), r0), r0_sq) + r2_axial = _div(_scale(_dot(r2, r0), r0), r0_sq) + r1_radial = _sub(r1, r1_axial) + r2_radial = _sub(r2, r2_axial) + r1_proj = _add(r1_axial, _div(_scale(epsilon, r1_radial), _norm(r1_radial))) + r2_proj = _add(r2_axial, _div(_scale(epsilon, r2_radial), _norm(r2_radial))) + r1Xr2_proj = _cross(r1_proj, r2_proj) + vel_ind_proj = _scale( gamma / (4 * np.pi) - * r1Xr2_proj - / (np.linalg.norm(r1Xr2_proj) ** 2) - * np.dot( + / (_norm(r1Xr2_proj) ** 2) + * _dot( r0, - r1_proj / np.linalg.norm(r1_proj) - r2_proj / np.linalg.norm(r2_proj), - ) + _sub(_div(r1_proj, _norm(r1_proj)), _div(r2_proj, _norm(r2_proj))), + ), + r1Xr2_proj, ) - return dist / epsilon * vel_ind_proj + return _scale(dist / epsilon, vel_ind_proj) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def _vel_trailing_vortex(XV1, XV2, XVP, gamma, Uinf): """Finite trailing-leg induced velocity (viscous core), as in Filament.velocity_3D_trailing_vortex.""" - r0 = XV2 - XV1 - r1 = XVP - XV1 - r2 = XVP - XV2 + r0 = _sub(XV2, XV1) + r1 = _sub(XVP, XV1) + r2 = _sub(XVP, XV2) - r_perp = r1 - np.dot(r1, r0) * r0 / (np.linalg.norm(r0) ** 2) - epsilon = np.sqrt(4 * _ALPHA0_OSEEN * _NU_AIR * np.linalg.norm(r_perp) / Uinf) + r0_sq = _norm(r0) ** 2 + r_perp = _sub(r1, _div(_scale(_dot(r1, r0), r0), r0_sq)) + epsilon = np.sqrt(4 * _ALPHA0_OSEEN * _NU_AIR * _norm(r_perp) / Uinf) - r1Xr0 = np.cross(r1, r0) - dist = np.linalg.norm(r1Xr0) / np.linalg.norm(r0) + r1Xr0 = _cross(r1, r0) + dist = _norm(r1Xr0) / _norm(r0) if dist > epsilon: - r1Xr2 = np.cross(r1, r2) - return ( + r1Xr2 = _cross(r1, r2) + return _scale( gamma / (4 * np.pi) - * r1Xr2 - / (np.linalg.norm(r1Xr2) ** 2) - * np.dot(r0, r1 / np.linalg.norm(r1) - r2 / np.linalg.norm(r2)) + / (_norm(r1Xr2) ** 2) + * _dot(r0, _sub(_div(r1, _norm(r1)), _div(r2, _norm(r2)))), + r1Xr2, ) elif dist < 1e-12 * epsilon: - return np.zeros(3) + return _ZERO3 else: - r1_radial = r1 - np.dot(r1, r0) * r0 / (np.linalg.norm(r0) ** 2) - r2_radial = r2 - np.dot(r2, r0) * r0 / (np.linalg.norm(r0) ** 2) - r1_proj = np.dot(r1, r0) * r0 / ( - np.linalg.norm(r0) ** 2 - ) + epsilon * r1_radial / np.linalg.norm(r1_radial) - r2_proj = np.dot(r2, r0) * r0 / ( - np.linalg.norm(r0) ** 2 - ) + epsilon * r2_radial / np.linalg.norm(r2_radial) - r1Xr2_proj = np.cross(r1_proj, r2_proj) - vel_ind_proj = ( + r1_axial = _div(_scale(_dot(r1, r0), r0), r0_sq) + r2_axial = _div(_scale(_dot(r2, r0), r0), r0_sq) + r1_radial = _sub(r1, r1_axial) + r2_radial = _sub(r2, r2_axial) + r1_proj = _add(r1_axial, _div(_scale(epsilon, r1_radial), _norm(r1_radial))) + r2_proj = _add(r2_axial, _div(_scale(epsilon, r2_radial), _norm(r2_radial))) + r1Xr2_proj = _cross(r1_proj, r2_proj) + vel_ind_proj = _scale( gamma / (4 * np.pi) - * r1Xr2_proj - / (np.linalg.norm(r1Xr2_proj) ** 2) - * np.dot( + / (_norm(r1Xr2_proj) ** 2) + * _dot( r0, - r1_proj / np.linalg.norm(r1_proj) - r2_proj / np.linalg.norm(r2_proj), - ) + _sub(_div(r1_proj, _norm(r1_proj)), _div(r2_proj, _norm(r2_proj))), + ), + r1Xr2_proj, ) - return dist / epsilon * vel_ind_proj + return _scale(dist / epsilon, vel_ind_proj) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def _vel_semiinfinite(XV1, Vf, XVP, gamma, Uinf, filament_direction): """Semi-infinite trailing-vortex induced velocity, as in Filament.velocity_3D_trailing_vortex_semiinfinite.""" GAMMA = -gamma * filament_direction - r1 = XVP - XV1 - r1XVf = np.cross(r1, Vf) - r_perp = r1 - np.dot(r1, Vf) * Vf + r1 = _sub(XVP, XV1) + r1XVf = _cross(r1, Vf) + r1_axial = _scale(_dot(r1, Vf), Vf) + r_perp = _sub(r1, r1_axial) - if np.linalg.norm(r_perp) < 1e-16: - return np.zeros(3) + if _norm(r_perp) < 1e-16: + return _ZERO3 - epsilon = np.sqrt(4 * _ALPHA0_OSEEN * _NU_AIR * np.linalg.norm(r_perp) / Uinf) + epsilon = np.sqrt(4 * _ALPHA0_OSEEN * _NU_AIR * _norm(r_perp) / Uinf) - dist = np.linalg.norm(r1XVf) / np.linalg.norm(Vf) + dist = _norm(r1XVf) / _norm(Vf) if dist > epsilon: K = ( GAMMA / 4 / np.pi - / np.linalg.norm(r1XVf) ** 2 - * (1 + np.dot(r1, Vf) / np.linalg.norm(r1)) + / _norm(r1XVf) ** 2 + * (1 + _dot(r1, Vf) / _norm(r1)) ) - return K * r1XVf + return _scale(K, r1XVf) elif dist < 1e-12 * epsilon: - return np.zeros(3) + return _ZERO3 else: - r1_radial = r1 - np.dot(r1, Vf) * Vf - r1_proj = np.dot(r1, Vf) * Vf + epsilon * r1_radial / np.linalg.norm(r1_radial) - r1XVf_proj = np.cross(r1_proj, Vf) + r1_radial = r_perp + r1_proj = _add(r1_axial, _div(_scale(epsilon, r1_radial), _norm(r1_radial))) + r1XVf_proj = _cross(r1_proj, Vf) K = ( GAMMA / 4 / np.pi - / np.linalg.norm(r1XVf_proj) ** 2 - * (1 + np.dot(r1_proj, Vf) / np.linalg.norm(r1_proj)) + / _norm(r1XVf_proj) ** 2 + * (1 + _dot(r1_proj, Vf) / _norm(r1_proj)) ) - return K * r1XVf_proj + return _scale(K, r1XVf_proj) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def _vel_bound_2D(XVP, bound_point_1, bound_point_2): """2D bound-vortex correction at the evaluation point, as in Panel.compute_velocity_induced_bound_2D (gamma = 1).""" - r3 = XVP - (bound_point_1 + bound_point_2) / 2 - r0 = bound_point_1 - bound_point_2 - r0_unit_X_r3 = np.cross(r0 / np.linalg.norm(r0), r3) - return (1 / (2 * np.pi)) * (r0_unit_X_r3 / (np.linalg.norm(r0_unit_X_r3) ** 2)) + r3 = _sub(XVP, _div(_add(bound_point_1, bound_point_2), 2.0)) + r0 = _sub(bound_point_1, bound_point_2) + r0_unit_X_r3 = _cross(_div(r0, _norm(r0)), r3) + return _scale(1 / (2 * np.pi), _div(r0_unit_X_r3, _norm(r0_unit_X_r3) ** 2)) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def assemble_AIC_matrices( eval_points, bound_point_1, @@ -212,43 +266,43 @@ def assemble_AIC_matrices( n = eval_points.shape[0] AIC = np.empty((3, n, n)) for icp in range(n): - ep = eval_points[icp] + ep = _row(eval_points, icp) for jring in range(n): + bp1 = _row(bound_point_1, jring) + bp2 = _row(bound_point_2, jring) + te1 = _row(TE_point_1, jring) + te2 = _row(TE_point_2, jring) if evaluation_point_on_bound: - velocity_induced = np.zeros(3) + velocity_induced = _ZERO3 else: velocity_induced = _vel_bound_vortex( - bound_point_2[jring], - bound_point_1[jring], - ep, - 1.0, - core_radius_fraction, + bp2, bp1, ep, 1.0, core_radius_fraction ) - wake_unit = wake_units[jring] + wake_unit = _row(wake_units, jring) wake_speed = wake_speeds[jring] - velocity_induced = velocity_induced + _vel_trailing_vortex( - bound_point_1[jring], TE_point_1[jring], ep, 1.0, wake_speed + velocity_induced = _add( + velocity_induced, _vel_trailing_vortex(bp1, te1, ep, 1.0, wake_speed) ) - velocity_induced = velocity_induced + _vel_trailing_vortex( - TE_point_2[jring], bound_point_2[jring], ep, 1.0, wake_speed + velocity_induced = _add( + velocity_induced, _vel_trailing_vortex(te2, bp2, ep, 1.0, wake_speed) ) - velocity_induced = velocity_induced + _vel_semiinfinite( - TE_point_1[jring], wake_unit, ep, 1.0, wake_speed, 1.0 + velocity_induced = _add( + velocity_induced, + _vel_semiinfinite(te1, wake_unit, ep, 1.0, wake_speed, 1.0), ) - velocity_induced = velocity_induced + _vel_semiinfinite( - TE_point_2[jring], wake_unit, ep, 1.0, wake_speed, -1.0 + velocity_induced = _add( + velocity_induced, + _vel_semiinfinite(te2, wake_unit, ep, 1.0, wake_speed, -1.0), ) if subtract_bound_2D and icp == jring: - velocity_induced = velocity_induced - _vel_bound_2D( - ep, bound_point_1[jring], bound_point_2[jring] - ) + velocity_induced = _sub(velocity_induced, _vel_bound_2D(ep, bp1, bp2)) AIC[0, icp, jring] = velocity_induced[0] AIC[1, icp, jring] = velocity_induced[1] AIC[2, icp, jring] = velocity_induced[2] return AIC -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def assemble_curved_bound_correction_AIC( eval_points, bound_point_1, @@ -315,50 +369,53 @@ def assemble_curved_bound_correction_AIC( for j in range(n): s_node[j + 1] = s_node[j] + L[j] for i in range(n): - ep = eval_points[i] - t = seg[i] / L[i] - c_i = chord_dir[i] + ep = _row(eval_points, i) + t = _div(_row(seg, i), L[i]) + c_i = _row(chord_dir, i) s_P = 0.5 * (s_node[i] + s_node[i + 1]) U = va_norm[i] for j in range(n): - ra = ep + sgn * (s_node[j] - s_P) * t - rb = ep + sgn * (s_node[j + 1] - s_P) * t + ra = _add(ep, _scale(sgn * (s_node[j] - s_P), t)) + rb = _add(ep, _scale(sgn * (s_node[j + 1] - s_P), t)) if fwd: rb2 = ra rb1 = rb else: rb2 = rb rb1 = ra - leg = (1.0 - ac_fraction) * chords[j] * c_i - rte1 = rb1 + leg - rte2 = rb2 + leg - bp1 = bound_point_1[j] - bp2 = bound_point_2[j] - te1 = TE_point_1[j] - te2 = TE_point_2[j] + leg = _scale((1.0 - ac_fraction) * chords[j], c_i) + rte1 = _add(rb1, leg) + rte2 = _add(rb2, leg) + bp1 = _row(bound_point_1, j) + bp2 = _row(bound_point_2, j) + te1 = _row(TE_point_1, j) + te2 = _row(TE_point_2, j) # unrolled legs (single filaments) minus actual legs (single filaments) - v = _vel_trailing_vortex(rb1, rte1, ep, 1.0, U) + _vel_trailing_vortex( - rte2, rb2, ep, 1.0, U + v = _add( + _vel_trailing_vortex(rb1, rte1, ep, 1.0, U), + _vel_trailing_vortex(rte2, rb2, ep, 1.0, U), ) - v = v - _vel_trailing_vortex(bp1, te1, ep, 1.0, U) - _vel_trailing_vortex( - te2, bp2, ep, 1.0, U + v = _sub( + _sub(v, _vel_trailing_vortex(bp1, te1, ep, 1.0, U)), + _vel_trailing_vortex(te2, bp2, ep, 1.0, U), ) # chordwise-spread actual system minus chordwise-spread unrolled system for q in range(n_strips): f = (strip_x[q] - ac_fraction) / (1.0 - ac_fraction) w = strip_w[q] - n1 = bp1 + f * (te1 - bp1) - n2 = bp2 + f * (te2 - bp2) - u1 = rb1 + f * (rte1 - rb1) - u2 = rb2 + f * (rte2 - rb2) - v = v + w * ( - _vel_bound_vortex(n2, n1, ep, 1.0, core_radius_fraction) - - _vel_bound_vortex(u2, u1, ep, 1.0, core_radius_fraction) - + _vel_trailing_vortex(n1, te1, ep, 1.0, U) - + _vel_trailing_vortex(te2, n2, ep, 1.0, U) - - _vel_trailing_vortex(u1, rte1, ep, 1.0, U) - - _vel_trailing_vortex(rte2, u2, ep, 1.0, U) + n1 = _add(bp1, _scale(f, _sub(te1, bp1))) + n2 = _add(bp2, _scale(f, _sub(te2, bp2))) + u1 = _add(rb1, _scale(f, _sub(rte1, rb1))) + u2 = _add(rb2, _scale(f, _sub(rte2, rb2))) + strip = _sub( + _vel_bound_vortex(n2, n1, ep, 1.0, core_radius_fraction), + _vel_bound_vortex(u2, u1, ep, 1.0, core_radius_fraction), ) + strip = _add(strip, _vel_trailing_vortex(n1, te1, ep, 1.0, U)) + strip = _add(strip, _vel_trailing_vortex(te2, n2, ep, 1.0, U)) + strip = _sub(strip, _vel_trailing_vortex(u1, rte1, ep, 1.0, U)) + strip = _sub(strip, _vel_trailing_vortex(rte2, u2, ep, 1.0, U)) + v = _add(v, _scale(w, strip)) A[0, i, j] = v[0] A[1, i, j] = v[1] A[2, i, j] = v[2] @@ -378,18 +435,18 @@ def thin_airfoil_strips(n_strips=24): return np.ascontiguousarray(x, dtype=float), np.ascontiguousarray(w / w.sum(), dtype=float) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def _on_segment_line(XV1, XV2, XVP): """True if XVP lies on the line through XV1 and XV2 (to round-off).""" - r0 = XV2 - XV1 - r1 = XVP - XV1 - len0 = np.linalg.norm(r0) + r0 = _sub(XV2, XV1) + r1 = _sub(XVP, XV1) + len0 = _norm(r0) if len0 == 0.0: return True - return np.linalg.norm(np.cross(r1, r0)) <= 1e-10 * len0 * max(np.linalg.norm(r1), len0) + return _norm(_cross(r1, r0)) <= 1e-10 * len0 * max(_norm(r1), len0) -@jit(nopython=True, cache=True) +@jit(nopython=True, cache=True, error_model="numpy") def induced_velocity_at_points( points, bound_point_1, @@ -413,31 +470,25 @@ def induced_velocity_at_points( n = bound_point_1.shape[0] out = np.zeros((m, 3)) for ip in range(m): - ep = points[ip] - vel = np.zeros(3) + ep = _row(points, ip) + vel = _ZERO3 for j in range(n): g = gamma[j] - wake_unit = wake_units[j] + bp1 = _row(bound_point_1, j) + bp2 = _row(bound_point_2, j) + te1 = _row(TE_point_1, j) + te2 = _row(TE_point_2, j) + wake_unit = _row(wake_units, j) wake_speed = wake_speeds[j] - vel = vel + _vel_bound_vortex( - bound_point_2[j], bound_point_1[j], ep, g, core_radius_fraction - ) + vel = _add(vel, _vel_bound_vortex(bp2, bp1, ep, g, core_radius_fraction)) # a straight filament induces nothing on its own line; skip the # legs the point lies on (the kernel would divide by zero there) - if not _on_segment_line(bound_point_1[j], TE_point_1[j], ep): - vel = vel + _vel_trailing_vortex( - bound_point_1[j], TE_point_1[j], ep, g, wake_speed - ) - if not _on_segment_line(TE_point_2[j], bound_point_2[j], ep): - vel = vel + _vel_trailing_vortex( - TE_point_2[j], bound_point_2[j], ep, g, wake_speed - ) - vel = vel + _vel_semiinfinite( - TE_point_1[j], wake_unit, ep, g, wake_speed, 1.0 - ) - vel = vel + _vel_semiinfinite( - TE_point_2[j], wake_unit, ep, g, wake_speed, -1.0 - ) + if not _on_segment_line(bp1, te1, ep): + vel = _add(vel, _vel_trailing_vortex(bp1, te1, ep, g, wake_speed)) + if not _on_segment_line(te2, bp2, ep): + vel = _add(vel, _vel_trailing_vortex(te2, bp2, ep, g, wake_speed)) + vel = _add(vel, _vel_semiinfinite(te1, wake_unit, ep, g, wake_speed, 1.0)) + vel = _add(vel, _vel_semiinfinite(te2, wake_unit, ep, g, wake_speed, -1.0)) out[ip, 0] = vel[0] out[ip, 1] = vel[1] out[ip, 2] = vel[2] diff --git a/tests/BodyAerodynamics/data/panel_kernel_reference.npz b/tests/BodyAerodynamics/data/panel_kernel_reference.npz new file mode 100644 index 0000000..897a218 Binary files /dev/null and b/tests/BodyAerodynamics/data/panel_kernel_reference.npz differ diff --git a/tests/BodyAerodynamics/panel_kernel_cases.py b/tests/BodyAerodynamics/panel_kernel_cases.py new file mode 100644 index 0000000..7a8431f --- /dev/null +++ b/tests/BodyAerodynamics/panel_kernel_cases.py @@ -0,0 +1,115 @@ +"""Geometry and solver cases for the force-law regression in +``test_panel_kernels.py``. + +``python tests/BodyAerodynamics/panel_kernel_cases.py OUT.npz`` writes every +numeric output of ``Solver.solve`` for every case. The committed reference +``data/panel_kernel_reference.npz`` was written with VSM at ``1201f6f``, the +last commit before the force law moved into ``VSM.core.panel_kernels``, so +the regression test pins the refactor to the original per-panel loop. +""" + +import sys +from pathlib import Path + +import numpy as np + +from VSM.core.BodyAerodynamics import BodyAerodynamics +from VSM.core.Solver import Solver +from VSM.core.WingGeometry import Wing + +REPO = Path(__file__).resolve().parents[2] +V3_FILE = REPO / "data/TUDELFT_V3_KITE/CAD_derived_geometry/aero_geometry_CAD_CFD_polars.yaml" +REFERENCE_FILE = Path(__file__).resolve().parent / "data" / "panel_kernel_reference.npz" + +_alpha = np.radians(np.linspace(-25, 25, 101)) +# Thin-airfoil lift with a small drag and moment, so all three magnitudes are exercised. +POLAR = np.column_stack( + (_alpha, 2 * np.pi * _alpha, 0.01 + 0.05 * _alpha**2, -0.1 - 0.02 * _alpha) +) + + +def _wing_body(leading_edges, chord_vector, n_panels): + wing = Wing(n_panels=n_panels, spanwise_panel_distribution="uniform") + for le in leading_edges: + wing.add_section(le, le + chord_vector, POLAR) + return BodyAerodynamics([wing]) + + +def straight_body(n_panels=16): + ys = np.linspace(-4.0, 4.0, n_panels + 1) + body = _wing_body([np.array([0.0, y, 0.0]) for y in ys], np.array([1.0, 0.0, 0.0]), n_panels) + body.va_initialize(20.0, 5.0, 0.0) + return body + + +def swept_body(n_panels=16): + """Symmetric 30 deg swept-back wing.""" + ys = np.linspace(-4.0, 4.0, n_panels + 1) + sweep = np.tan(np.radians(30.0)) + body = _wing_body( + [np.array([sweep * abs(y), y, 0.0]) for y in ys], np.array([1.0, 0.0, 0.0]), n_panels + ) + body.va_initialize(15.0, 6.0, 2.0) + return body + + +def arc_body(n_panels=16): + """Circular-arc wing (kite-like anhedral), constant chord along x.""" + th = np.radians(np.linspace(-55.0, 55.0, n_panels + 1)) + les = [np.array([0.0, 4.0 * np.sin(t), 4.0 * np.cos(t)]) for t in th] + body = _wing_body(les, np.array([1.0, 0.0, 0.0]), n_panels) + body.va_initialize(10.0, 8.0, 0.0) + return body + + +def v3_body(n_panels=20, yaw_rate=0.0): + body = BodyAerodynamics.instantiate(n_panels=n_panels, file_path=V3_FILE) + body.va_initialize(10.0, 8.0, 3.0, np.array([0.0, 0.0, yaw_rate])) + return body + + +# name -> (body factory, solver kwargs) +CASES = { + "straight": (straight_body, {}), + "swept": (swept_body, {}), + "arc": (arc_body, {}), + "v3": (v3_body, {}), + "v3_viscous_noAT": ( + v3_body, + {"is_with_viscous_drag_correction": True, "is_with_attached_trailed_vortex_force": False}, + ), + "v3_yaw_rate": (lambda: v3_body(yaw_rate=0.4), {"reference_point": [0.1, 0.0, 0.5]}), +} + + +def numeric_outputs(results: dict) -> dict: + """Every numeric entry of a ``solve`` result, as float arrays.""" + out = {} + for key, value in results.items(): + try: + array = np.asarray(value, dtype=float) + except (TypeError, ValueError): + continue + out[key] = array + return out + + +def run_case(name: str, is_aoa_corrected: bool) -> dict: + factory, kwargs = CASES[name] + solver = Solver(allowed_error=1e-10, is_aoa_corrected=is_aoa_corrected, **kwargs) + results = solver.solve(factory()) + assert results["gamma_converged"] + return numeric_outputs(results) + + +def all_outputs() -> dict: + flat = {} + for name in CASES: + for corrected in (True, False): + for key, value in run_case(name, corrected).items(): + flat[f"{name}|{corrected}|{key}"] = value + return flat + + +if __name__ == "__main__": + np.savez(sys.argv[1], **all_outputs()) diff --git a/tests/BodyAerodynamics/test_panel_kernels.py b/tests/BodyAerodynamics/test_panel_kernels.py new file mode 100644 index 0000000..cc20943 --- /dev/null +++ b/tests/BodyAerodynamics/test_panel_kernels.py @@ -0,0 +1,303 @@ +"""The section force law as ``xp`` kernels (``VSM.core.panel_kernels``). + +Contract: + +* moving the force law into the kernels changed no output: every numeric + result of ``Solver.solve`` equals the one of VSM ``1201f6f`` (the commit + before the move) to round-off, for ``is_aoa_corrected`` True and False, on + a straight, a swept, a curved (arc) wing and the LEI V3 kite; +* the public quarter-chord induction matrices are the ones the force + directions are computed from; +* every kernel gives the same numbers with ``xp=casadi`` as with + ``xp=numpy``; +* a CasADi graph with the induction matrices as parameters reproduces the + per-panel forces of ``compute_results``. +""" + +import sys +from pathlib import Path + +import numpy as np +import pytest + +sys.path.insert(0, str(Path(__file__).resolve().parent)) + +from panel_kernel_cases import CASES, REFERENCE_FILE, run_case, v3_body, arc_body # noqa: E402 +from VSM.core import panel_kernels # noqa: E402 +from VSM.core.Solver import Solver # noqa: E402 +from VSM.core.utils import thin_airfoil_strips, assemble_curved_bound_correction_AIC # noqa: E402 + +RTOL = 1e-12 + +_FORCE_KEYS = ("Fx", "Fy", "Fz") +_MOMENT_KEYS = ("Mx", "My", "Mz") + + +# ----------------------------------------------------------- no-change regression +@pytest.fixture(scope="module") +def reference(): + return np.load(REFERENCE_FILE) + + +@pytest.mark.parametrize("is_aoa_corrected", [True, False]) +@pytest.mark.parametrize("case", list(CASES)) +def test_solve_outputs_match_reference_commit(reference, case, is_aoa_corrected): + """Every numeric output equals VSM 1201f6f to round-off. Entries that + vanish by symmetry (side force of a symmetric wing, ...) are compared + against the case's force or moment scale instead of their own.""" + outputs = run_case(case, is_aoa_corrected) + prefix = f"{case}|{is_aoa_corrected}|" + expected_keys = {k[len(prefix):] for k in reference.files if k.startswith(prefix)} + assert set(outputs) == expected_keys + family_scale = { + "force": max(abs(float(reference[prefix + k])) for k in _FORCE_KEYS), + "moment": max(abs(float(reference[prefix + k])) for k in _MOMENT_KEYS), + "other": max(abs(float(reference[prefix + k])) for k in ("cl", "cd")), + } + for key, value in outputs.items(): + ref = reference[prefix + key] + assert value.shape == ref.shape, key + if not ref.size: + continue + if key[0] == "F" or key in ("lift", "drag", "side", "drag_induced_trefftz", "bridle_line_forces"): + family = "force" + elif key[0] == "M": + family = "moment" + else: + family = "other" + scale = max(np.max(np.abs(ref)), 1e-3 * family_scale[family]) + np.testing.assert_allclose(value, ref, rtol=0.0, atol=RTOL * scale, err_msg=key) + + +# ------------------------------------------------------ quarter-chord induction +@pytest.mark.parametrize("factory", [arc_body, v3_body]) +def test_quarter_chord_matrices_are_llt_plus_curved_correction(factory): + """The public matrices equal the trailed system on the quarter chord plus + the curved bound-vortex correction, assembled independently here.""" + body = factory() + va = np.array([p.va for p in body.panels]) + va_norm = np.linalg.norm(va, axis=1) + va_unit = va / va_norm[:, None] + core = 0.05 + public = body.compute_quarter_chord_AIC_matrices(core, va_norm, va_unit) + + llt = body.compute_AIC_matrices("LLT", core, va_norm, va_unit) + rows = lambda name: np.ascontiguousarray([getattr(p, name) for p in body.panels], dtype=float) # noqa: E731 + strip_x, strip_w = thin_airfoil_strips() + correction = assemble_curved_bound_correction_AIC( + rows("aerodynamic_center"), rows("bound_point_1"), rows("bound_point_2"), + rows("TE_point_1"), rows("TE_point_2"), rows("y_airf"), + np.array([p.chord for p in body.panels]), va_norm, + float(body._aerodynamic_center_location), core, strip_x, strip_w, + ) + for k in range(3): + assert public[k].shape == (body.n_panels, body.n_panels) + np.testing.assert_allclose(public[k], llt[k] + correction[k], rtol=1e-13, atol=1e-15) + if factory is arc_body: # curved: the correction is not negligible + assert np.max(np.abs(correction)) > 1e-3 * np.max(np.abs(llt)) + + +def test_quarter_chord_direction_angle_uses_the_public_matrices(): + body = v3_body() + solver = Solver(allowed_error=1e-10) + results = solver.solve(body) + va_norm = np.linalg.norm(solver.va_array, axis=1) + AIC_qc = body.compute_quarter_chord_AIC_matrices( + solver.core_radius_fraction, va_norm, solver.va_array / va_norm[:, None] + ) + gamma = results["gamma_distribution"] + v_qc = solver.va_array + np.column_stack([A @ gamma for A in AIC_qc]) + alpha_qc = np.arctan2( + np.sum(solver.x_airf_array * v_qc, axis=1), np.sum(solver.y_airf_array * v_qc, axis=1) + ) + np.testing.assert_allclose(np.ravel(results["alpha_at_ac"]), alpha_qc, rtol=0, atol=1e-14) + + +# --------------------------------------------------------- casadi == numpy +@pytest.fixture(scope="module") +def solved_v3(): + """A converged V3 solve and the arrays a force-law caller needs.""" + body = v3_body() + solver = Solver( + allowed_error=1e-10, + is_aoa_corrected=True, + is_with_attached_trailed_vortex_force=True, + ) + results = solver.solve(body) + va_norm = np.linalg.norm(solver.va_array, axis=1) + AIC_qc = body.compute_quarter_chord_AIC_matrices( + solver.core_radius_fraction, va_norm, solver.va_array / va_norm[:, None] + ) + return body, solver, results, AIC_qc + + +def _numpy_inputs(solver, results, AIC_qc): + gamma = np.asarray(results["gamma_distribution"], dtype=float) + aic = (solver.AIC_x, solver.AIC_y, solver.AIC_z) + v_rel = panel_kernels.quarter_chord_velocity(solver.va_array, gamma, *aic) + alpha, umag = panel_kernels.relative_flow( + v_rel, solver.x_airf_array, solver.y_airf_array, solver.z_airf_array + ) + panels = solver.panels + return { + "gamma": gamma, + "v_rel": v_rel, + "alpha": alpha, + "umag": umag, + "cl": np.array([p.compute_cl(a) for p, a in zip(panels, alpha)]), + "cd": np.array([p.compute_cd_cm(a)[0] for p, a in zip(panels, alpha)]), + "cm": np.array([p.compute_cd_cm(a)[1] for p, a in zip(panels, alpha)]), + "AIC_qc": AIC_qc, + } + + +def test_kernels_casadi_equals_numpy(solved_v3): + ca = pytest.importorskip("casadi") + body, solver, results, AIC_qc = solved_v3 + d = _numpy_inputs(solver, results, AIC_qc) + x, y, z = solver.x_airf_array, solver.y_airf_array, solver.z_airf_array + chord, width, rho = solver.chord_array, solver.width_array, solver.rho + + col = lambda v: ca.DM(np.asarray(v, dtype=float).reshape(-1, 1)) # noqa: E731 + mat = lambda a: ca.DM(np.asarray(a, dtype=float)) # noqa: E731 + + def same(ca_value, np_value): + ca_array = np.array(ca_value, dtype=float).reshape(np.shape(np_value)) + scale = max(np.max(np.abs(np_value)), 1e-300) + np.testing.assert_allclose(ca_array, np_value, rtol=0, atol=RTOL * scale) + + # flow + same(panel_kernels.relative_flow(mat(d["v_rel"]), mat(x), mat(y), mat(z), xp=ca)[0], d["alpha"]) + same(panel_kernels.relative_flow(mat(d["v_rel"]), mat(x), mat(y), mat(z), xp=ca)[1], d["umag"]) + v_qc_np = panel_kernels.quarter_chord_velocity(solver.va_array, d["gamma"], *AIC_qc) + v_qc_ca = panel_kernels.quarter_chord_velocity( + mat(solver.va_array), col(d["gamma"]), *[mat(A) for A in AIC_qc], xp=ca + ) + same(v_qc_ca, v_qc_np) + alpha_dir_np = panel_kernels.direction_angle(v_qc_np, x, y) + same(panel_kernels.direction_angle(v_qc_ca, mat(x), mat(y), xp=ca), alpha_dir_np) + + # force law, both direction laws + for alpha_dir in (None, alpha_dir_np): + np_out = panel_kernels.panel_forces( + d["alpha"], d["umag"], d["cl"], d["cd"], d["cm"], chord, width, x, y, z, rho, + alpha_dir=alpha_dir, + ) + ca_out = panel_kernels.panel_forces( + col(d["alpha"]), col(d["umag"]), col(d["cl"]), col(d["cd"]), col(d["cm"]), + col(chord), col(width), mat(x), mat(y), mat(z), rho, + alpha_dir=None if alpha_dir is None else col(alpha_dir), xp=ca, + ) + for ca_value, np_value in zip(ca_out, np_out): + same(ca_value, np_value) + np_vec = panel_kernels.section_force_vectors( + d["alpha"], d["umag"], d["cl"], d["cd"], d["cm"], chord, x, y, z, rho, + alpha_dir=alpha_dir, + ) + ca_vec = panel_kernels.section_force_vectors( + col(d["alpha"]), col(d["umag"]), col(d["cl"]), col(d["cd"]), col(d["cm"]), + col(chord), mat(x), mat(y), mat(z), rho, + alpha_dir=None if alpha_dir is None else col(alpha_dir), xp=ca, + ) + for ca_value, np_value in zip(ca_vec, np_vec): + same(ca_value, np_value) + + # attached trailed vortex force + legs = np.random.default_rng(0).normal(size=(body.n_panels + 1, 3)) + v_ind = np.random.default_rng(1).normal(size=(body.n_panels + 1, 3)) + same( + panel_kernels.attached_trailed_forces( + col(d["gamma"]), mat(solver.va_array), mat(v_ind), legs, rho, xp=ca + ), + panel_kernels.attached_trailed_forces(d["gamma"], solver.va_array, v_ind, legs, rho), + ) + + +# ------------------------------------------- casadi graph == compute_results +def _boundary_shares(n): + """(n, n+1) split of boundary loads onto panels, compute_results' rule.""" + shares = np.zeros((n, n + 1)) + shares[0, 0] = 1.0 + shares[n - 1, n] = 1.0 + for j in range(1, n): + shares[j - 1, j] = shares[j, j] = 0.5 + return shares + + +def _compute_results_at(body, solver, gamma): + """compute_results for this ``gamma`` with ``alpha``/``umag`` evaluated at + the same ``gamma``. (``solve`` passes the loop's last alpha, which lags the + returned gamma by one iteration, i.e. agrees with it only to the + convergence tolerance.)""" + alpha, umag, _ = solver.compute_aerodynamic_quantities(gamma) + va_norm = np.linalg.norm(solver.va_array, axis=1) + return body.compute_results( + gamma, solver.rho, solver.aerodynamic_model_type, solver.core_radius_fraction, + solver.mu, alpha, umag, solver.chord_array, solver.x_airf_array, + solver.y_airf_array, solver.z_airf_array, solver.va_array, va_norm, + solver.va_array / va_norm[:, None], solver.panels, False, + solver.is_with_viscous_drag_correction, solver.reference_point, + solver.is_aoa_corrected, + relative_velocity_array=solver.compute_relative_velocity(gamma), + is_with_attached_trailed_vortex_force=solver.is_with_attached_trailed_vortex_force, + ) + + +@pytest.mark.parametrize("is_aoa_corrected", [True, False]) +def test_casadi_graph_reproduces_compute_results_forces(is_aoa_corrected): + """Graph: gamma, va and the 3/4-chord and quarter-chord matrices are MX + parameters; ``v_rel`` and ``v_qc`` come from ``quarter_chord_velocity``, + the section forces from ``Solver.casadi_panel_force_function``. Evaluated + at a converged gamma it reproduces compute_results' panel forces and + moments (the attached-trailed part removed, as it is not in the + function).""" + ca = pytest.importorskip("casadi") + body = v3_body() + solver = Solver(allowed_error=1e-10, is_aoa_corrected=is_aoa_corrected) + gamma_converged = solver.solve(body)["gamma_distribution"] + results = _compute_results_at(body, solver, gamma_converged) + n = solver.n_panels + va_norm = np.linalg.norm(solver.va_array, axis=1) + AIC_qc = body.compute_quarter_chord_AIC_matrices( + solver.core_radius_fraction, va_norm, solver.va_array / va_norm[:, None] + ) + + gamma = ca.MX.sym("gamma", n) + va = ca.MX.sym("va", n, 3) + A_cp = [ca.MX.sym(f"A_cp_{k}", n, n) for k in range(3)] + A_qc = [ca.MX.sym(f"A_qc_{k}", n, n) for k in range(3)] + v_rel = panel_kernels.quarter_chord_velocity(va, gamma, *A_cp, xp=ca) + v_qc = panel_kernels.quarter_chord_velocity(va, gamma, *A_qc, xp=ca) if is_aoa_corrected else v_rel + forces_fn = solver.casadi_panel_force_function() + force, moment_local = forces_fn( + v_rel, v_qc, solver.x_airf_array, solver.y_airf_array, solver.z_airf_array, + solver.chord_array, solver.width_array, solver.rho, + )[:2] + graph = ca.Function("forces", [gamma, va, *A_cp, *A_qc], [force, moment_local]) + + force_num, moment_num = graph( + gamma_converged, solver.va_array, + solver.AIC_x, solver.AIC_y, solver.AIC_z, *AIC_qc, + ) + force_num = np.array(force_num) + moment_num = np.array(moment_num) + + shares = _boundary_shares(n) + f_at = np.array(results["F_attached_trailed_distribution"]) + assert f_at.shape == (n + 1, 3) # AT force on (the default) + expected_force = np.array(results["F_distribution"]) - shares @ f_at + reference_point = np.zeros(3) + ac = np.array([p.aerodynamic_center for p in solver.panels]) + at_points = body.compute_attached_trailed_vortex_forces( + gamma_converged, solver.rho, solver.core_radius_fraction, solver.va_array + )[1] + expected_moment = ( + np.array(results["M_distribution"]) + - shares @ np.cross(at_points - reference_point, f_at) + - np.cross(ac - reference_point, expected_force) + ) + f_scale = np.max(np.abs(expected_force)) + np.testing.assert_allclose(force_num, expected_force, rtol=0, atol=RTOL * f_scale) + np.testing.assert_allclose( + moment_num, expected_moment, rtol=0, atol=RTOL * np.max(np.abs(expected_moment)) + )