diff --git a/src/VSM/core/Solver.py b/src/VSM/core/Solver.py index 96d3766..bb03468 100644 --- a/src/VSM/core/Solver.py +++ b/src/VSM/core/Solver.py @@ -6,6 +6,26 @@ from VSM.core import panel_kernels + +def _table_symbols(n_panels: int, alpha_tables) -> dict: + """Symbolic ``(n_panels, n_alpha)`` polar tables for ``table_inputs``.""" + import casadi as ca + + n_alpha = int(np.asarray(alpha_tables[0]).size) + return {name: ca.MX.sym(f"{name}_table", n_panels, n_alpha) for name in ("cl", "cd", "cm")} + + +def _with_table_inputs(inputs: list, input_names: list, table_syms): + """Append the table inputs (if any) and default them to NaN, so a call + that omits them (CasADi fills omitted inputs with 0 by default) gives NaN + rather than a silently zero polar.""" + if table_syms is None: + return inputs, input_names, {} + inputs = inputs + [table_syms[k] for k in ("cl", "cd", "cm")] + input_names = input_names + ["cl_table", "cd_table", "cm_table"] + defaults = [0.0] * (len(inputs) - 3) + [float("nan")] * 3 + return inputs, input_names, {"default_in": defaults} + class Solver: """Solver for aerodynamic circulation distribution and force computation. @@ -203,6 +223,7 @@ def __init__( self._newton_reject_ratio = 0.5 self._casadi_newton_cache: dict = {} self._casadi_newton_loop_cache: dict = {} + self._casadi_table_input_cache: dict = {} # polar tables as inputs #: 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. @@ -1030,12 +1051,91 @@ def _casadi_newton_key(self, alpha_tables, cl_tables, cd_tables, cm_tables) -> t digest.hexdigest(), ) - def _casadi_polar_lookups(self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None): + def _casadi_table_input_key(self, alpha_tables) -> tuple: + """Cache key of the functions that take the polar tables as INPUTS + (``table_inputs=True``): the graph depends on the panel count, the + shared alpha grid, the interpolation and the regularization settings, + not on the table values, so one function serves every shape whose + blended polars change but keep the grid.""" + self._check_table_inputs_supported(alpha_tables) + grid = np.ascontiguousarray(alpha_tables[0], dtype=float) + return ( + "table_inputs", + self.n_panels, + bool(self.is_with_artificial_viscosity), + float(self.artificial_viscosity_factor), + self.polar_interpolation, + hashlib.blake2b(grid.tobytes(), digest_size=16).hexdigest(), + ) + + def _check_table_inputs_supported(self, alpha_tables) -> None: + if self.polar_interpolation != "linear": + # A parametric cubic spline refits its coefficients on every call + # (~0.1 s per evaluation), which defeats the purpose. + raise ValueError( + "Polar tables as CasADi inputs need polar_interpolation='linear' " + f"(got {self.polar_interpolation!r})." + ) + if not all( + table.shape == alpha_tables[0].shape and np.array_equal(table, alpha_tables[0]) + for table in alpha_tables[1:] + ): + raise ValueError( + "Polar tables as CasADi inputs need every panel polar on one " + "shared alpha grid." + ) + + def casadi_polar_table_inputs(self) -> dict: + """The current panels' polar tables in the layout of the ``table_inputs`` + CasADi functions (:meth:`_casadi_newton_function` and + :meth:`casadi_panel_force_function` with ``table_inputs=True``): + ``{"cl_table", "cd_table", "cm_table"}``, each ``(n_panels, n_alpha)``, + row ``i`` = panel ``i``'s coefficients on the shared alpha grid + (``panel.panel_polar_data[:, 0]``, radians). + + Pass them on every call, next to the geometry: a deformed shape blends + its section polars with new weights, so the values change while the + grid, and with it the function, stays the same. + """ + alpha_tables, cl_tables, cd_tables, cm_tables = self._polar_tables() + self._check_table_inputs_supported(alpha_tables) + return { + "cl_table": np.vstack(cl_tables), + "cd_table": np.vstack(cd_tables), + "cm_table": np.vstack(cm_tables), + } + + def casadi_stall_angle_inputs(self) -> dict: + """``{"stall_angles", "stall_angles_neg"}`` [rad] of the current panel + polars, the gate inputs of the residual function as + :meth:`gamma_loop_casadi_newton` passes them (plus / minus infinity + without artificial viscosity). They are derived from the tables, so a + caller passing new tables passes these new values with them.""" + if self.is_with_artificial_viscosity: + return { + "stall_angles": self._panel_stall_angles(), + "stall_angles_neg": self._panel_negative_stall_angles(), + } + return { + "stall_angles": np.full(self.n_panels, np.inf), + "stall_angles_neg": np.full(self.n_panels, -np.inf), + } + + def _casadi_polar_lookups( + self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None, table_symbols=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. + + ``table_symbols`` (``{"cl": MX, "cd": MX, "cm": MX}``, each + ``(n, n_alpha)``) makes the table VALUES symbolic inputs instead of + constants: one 1-D parametric linear interpolant on the shared grid, + mapped over the panels. It reproduces the constant table exactly + (same nodes, same clamping); only the shared-grid linear case is + supported (see :meth:`_check_table_inputs_supported`). """ import casadi as ca @@ -1059,6 +1159,33 @@ def _casadi_polar_lookups(self, alpha_tables, cl_tables, cd_tables=None, cm_tabl kind = self.polar_interpolation smooth = kind == "bspline" + if table_symbols is not None: + self._check_table_inputs_supported(alpha_tables) + grid = np.asarray(alpha_tables[0], dtype=float) + lo, hi = float(grid[0]), float(grid[-1]) + row_interpolant = ca.interpolant("polar_row", "linear", [grid], 1) + alpha_row = ca.MX.sym("alpha_row") + values_row = ca.MX.sym("values_row", grid.size) + # One panel's lookup; .map(n) evaluates it column-wise over panels. + row_lookup = ca.Function( + "polar_row_lookup", + [alpha_row, values_row], + [row_interpolant(ca.fmin(ca.fmax(alpha_row, lo), hi), values_row)], + ).map(n) + + def make_parametric_lookup(table): + def lookup(alpha_vec): + return row_lookup(alpha_vec.T, table.T).T + + return lookup + + return ( + make_parametric_lookup(table_symbols["cl"]), + make_parametric_lookup(table_symbols["cd"]), + make_parametric_lookup(table_symbols["cm"]), + smooth, + ) + def make_lookup(name, value_tables): """Clamped lookup ``alpha_vec (n) -> values (n)`` of one polar column.""" if shared_grid: @@ -1104,7 +1231,9 @@ def lookup(alpha_vec): 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): + def _build_casadi_newton_function( + self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None, table_inputs=False + ): """Build the CasADi function of the PER-PANEL section physics that :meth:`gamma_loop_casadi_newton` assembles its residual and exact Jacobian from. @@ -1154,12 +1283,18 @@ def _build_casadi_newton_function(self, alpha_tables, cl_tables, cd_tables=None, 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. + + ``table_inputs=True`` appends the inputs ``cl_table``, ``cd_table``, + ``cm_table`` (each ``(n, n_alpha)``, :meth:`casadi_polar_table_inputs`) + and takes the table values from them instead of embedding them; see + :meth:`_casadi_newton_function`. """ import casadi as ca n = self.n_panels + table_syms = _table_symbols(n, alpha_tables) if table_inputs else None cl_of, cd_of, cm_of, smooth = self._casadi_polar_lookups( - alpha_tables, cl_tables, cd_tables, cm_tables + alpha_tables, cl_tables, cd_tables, cm_tables, table_symbols=table_syms ) v_rel = ca.MX.sym("v_rel", n, 3) @@ -1244,9 +1379,12 @@ def averaged(lookup): seeds = [ca.DM(np.eye(3)[k][None, :].repeat(n, axis=0)) for k in range(3)] dh_dvrel = ca.horzcat(*[ca.jtimes(h_for_jacobian, v_rel, seed) for seed in seeds]) + inputs = [v_rel, x_airf, y_airf, z_airf, chord, width, stall_angles, stall_angles_neg, l_gamma] + input_names = ["v_rel", "x_airf", "y_airf", "z_airf", "chord", "width", "stall_angles", "stall_angles_neg", "l_gamma"] + inputs, input_names, options = _with_table_inputs(inputs, input_names, table_syms) return ca.Function( - "vsm_panel_residual", - [v_rel, x_airf, y_airf, z_airf, chord, width, stall_angles, stall_angles_neg, l_gamma], + "vsm_panel_residual_tables" if table_inputs else "vsm_panel_residual", + inputs, [ h, dh_dvrel, @@ -1262,11 +1400,29 @@ def averaged(lookup): averaged(cd_of), averaged(cm_of), ], - ["v_rel", "x_airf", "y_airf", "z_airf", "chord", "width", "stall_angles", "stall_angles_neg", "l_gamma"], + input_names, ["h", "dh_dvrel", "mu", "alpha", "umag", "gate", "h_avg", "cl", "cd", "cm", "cl_avg", "cd_avg", "cm_avg"], + options, ) - def _casadi_newton_function(self): + def _casadi_newton_function(self, table_inputs: bool = False): + """The per-panel residual function (:meth:`_build_casadi_newton_function`) + of the current panels, cached. + + ``table_inputs=False`` (default, what :meth:`gamma_loop_casadi_newton` + uses): the polar tables are constants of the graph, cached per table + set. ``table_inputs=True``: the tables are the extra inputs + ``cl_table, cd_table, cm_table``, cached per alpha grid, so a caller + that re-trims deformed shapes builds the function once and passes + :meth:`casadi_polar_table_inputs` (and :meth:`casadi_stall_angle_inputs`) + per shape. Same outputs to round-off; needs ``polar_interpolation= + "linear"`` and one shared alpha grid. Omitted table inputs default to + NaN, so a call that forgets them fails loudly instead of using zeros. + """ + if table_inputs: + return self._cached_table_input_function( + "newton", self._build_casadi_newton_function + ) alpha_tables, cl_tables, cd_tables, cm_tables = self._polar_tables() key = self._casadi_newton_key(alpha_tables, cl_tables, cd_tables, cm_tables) fn = self._casadi_newton_cache.get(key) @@ -1282,7 +1438,7 @@ def _casadi_newton_function(self): return fn def _build_casadi_panel_force_function( - self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None + self, alpha_tables, cl_tables, cd_tables=None, cm_tables=None, table_inputs=False ): """CasADi function of the per-panel section forces, the force law of ``compute_results`` (``VSM.core.panel_kernels.panel_forces``) on the @@ -1291,8 +1447,9 @@ def _build_casadi_panel_force_function( import casadi as ca n = self.n_panels + table_syms = _table_symbols(n, alpha_tables) if table_inputs else None cl_of, cd_of, cm_of, _ = self._casadi_polar_lookups( - alpha_tables, cl_tables, cd_tables, cm_tables + alpha_tables, cl_tables, cd_tables, cm_tables, table_symbols=table_syms ) v_rel = ca.MX.sym("v_rel", n, 3) v_qc = ca.MX.sym("v_qc", n, 3) @@ -1310,15 +1467,21 @@ def _build_casadi_panel_force_function( 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", + inputs, input_names, options = _with_table_inputs( [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"], + table_syms, + ) + return ca.Function( + "vsm_panel_forces_tables" if table_inputs else "vsm_panel_forces", + inputs, + [force, moment_local, alpha, alpha_dir, umag, cl, cd, cm], + input_names, ["force", "moment_local", "alpha", "alpha_dir", "umag", "cl", "cd", "cm"], + options, ) - def casadi_panel_force_function(self): + def casadi_panel_force_function(self, table_inputs: bool = False): """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. @@ -1343,9 +1506,16 @@ def casadi_panel_force_function(self): 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. + + ``table_inputs=True`` appends ``cl_table, cd_table, cm_table`` as + inputs (see :meth:`_casadi_newton_function`). """ import casadi # noqa: F401 (clear error when the optional dep is missing) + if table_inputs: + return self._cached_table_input_function( + "panel_forces", self._build_casadi_panel_force_function + ) 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 @@ -1360,6 +1530,20 @@ def casadi_panel_force_function(self): self._casadi_newton_cache[key] = fn return fn + def _cached_table_input_function(self, kind: str, builder): + """Function of ``kind`` with the polar tables as inputs, built once + per (panel count, alpha grid, settings) and then reused for any + table values.""" + alpha_tables, cl_tables, cd_tables, cm_tables = self._polar_tables() + key = (kind,) + self._casadi_table_input_key(alpha_tables) + fn = self._casadi_table_input_cache.get(key) + if fn is None: + fn = builder(alpha_tables, cl_tables, cd_tables, cm_tables, table_inputs=True) + if len(self._casadi_table_input_cache) >= 32: + self._casadi_table_input_cache.pop(next(iter(self._casadi_table_input_cache))) + self._casadi_table_input_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 diff --git a/tests/Solver/test_casadi_polar_table_inputs.py b/tests/Solver/test_casadi_polar_table_inputs.py new file mode 100644 index 0000000..34ed806 --- /dev/null +++ b/tests/Solver/test_casadi_polar_table_inputs.py @@ -0,0 +1,168 @@ +"""Polar tables as inputs of the CasADi panel functions (``table_inputs=True``). + +Contract: with the tables passed as inputs, the residual and panel-force +functions give the same outputs as the functions with the tables embedded; +one function, built once, serves a second shape whose blended polars differ +(no rebuild); omitted tables give NaN, not a silently zero polar; and the +function still expands to SX. +""" + +import sys +from pathlib import Path + +import numpy as np +import pytest + +ca = pytest.importorskip("casadi") + +REPO = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(REPO / "tests" / "BodyAerodynamics")) + +from panel_kernel_cases import V3_FILE # noqa: E402 +from VSM.core.BodyAerodynamics import BodyAerodynamics # noqa: E402 +from VSM.core.Solver import Solver # noqa: E402 + +RTOL = 1e-12 + + +def _v3(distribution): + body = BodyAerodynamics.instantiate( + n_panels=20, file_path=V3_FILE, spanwise_panel_distribution=distribution + ) + body.va_initialize(10.0, 8.0, 3.0) + return body + + +def _solved(body, **solver_kwargs): + solver = Solver(gamma_loop_type="casadi_newton", allowed_error=1e-10, **solver_kwargs) + results = solver.solve(body) + return solver, np.asarray(results["gamma_distribution"], dtype=float) + + +def _newton_inputs(solver, gamma): + v_rel = solver.compute_relative_velocity(gamma) + n = solver.n_panels + laplacian = ( + solver._build_spanwise_laplacian() + if solver.is_with_artificial_viscosity + else np.zeros((n, n)) + ) + return { + "v_rel": v_rel, + "x_airf": solver.x_airf_array, + "y_airf": solver.y_airf_array, + "z_airf": solver.z_airf_array, + "chord": solver.chord_array, + "width": solver.width_array, + "l_gamma": laplacian @ gamma, + **solver.casadi_stall_angle_inputs(), + } + + +def _force_inputs(solver, gamma): + v_rel = solver.compute_relative_velocity(gamma) + return { + "v_rel": v_rel, + "v_qc": v_rel, + "x_airf": solver.x_airf_array, + "y_airf": solver.y_airf_array, + "z_airf": solver.z_airf_array, + "chord": solver.chord_array, + "width": solver.width_array, + "rho": solver.rho, + } + + +def _assert_outputs_equal(out_tables, out_const): + assert set(out_tables) == set(out_const) + for key, value in out_const.items(): + expected = np.array(value, dtype=float) + actual = np.array(out_tables[key], dtype=float) + scale = max(np.max(np.abs(expected)), 1e-300) + np.testing.assert_allclose(actual, expected, rtol=0, atol=RTOL * scale, err_msg=key) + + +@pytest.mark.parametrize("viscosity", [False, True]) +def test_table_inputs_reproduce_embedded_tables(viscosity): + """Residual function (all 13 outputs, incl. the corner-averaged ones and + the viscosity lift slope) and panel-force function, against the + embedded-table functions of the same shape.""" + solver, gamma = _solved(_v3("uniform"), is_with_artificial_viscosity=viscosity) + tables = solver.casadi_polar_table_inputs() + assert tables["cl_table"].shape == (solver.n_panels, solver.panels[0].panel_polar_data.shape[0]) + + args = _newton_inputs(solver, gamma) + if viscosity: # put a panel past stall so the gate and mu are exercised + args["stall_angles"] = args["stall_angles"] - 0.3 + _assert_outputs_equal( + solver._casadi_newton_function(table_inputs=True)(**args, **tables), + solver._casadi_newton_function()(**args), + ) + f_args = _force_inputs(solver, gamma) + _assert_outputs_equal( + solver.casadi_panel_force_function(table_inputs=True)(**f_args, **tables), + solver.casadi_panel_force_function()(**f_args), + ) + + +def test_one_function_serves_a_second_shape_without_rebuild(): + """Shapes A and B have the same panel count and alpha grid but differently + blended polars. The table-input function built for A, given B's tables and + geometry, equals a fresh embedded-table build for B, and no new + table-input function is built.""" + solver_a, gamma_a = _solved(_v3("uniform")) + fn = solver_a._casadi_newton_function(table_inputs=True) + force_fn = solver_a.casadi_panel_force_function(table_inputs=True) + cache_size = len(solver_a._casadi_table_input_cache) + + body_b = _v3("cosine") + tables_a = solver_a.casadi_polar_table_inputs() + solver_a.solve(body_b) # same solver object, now on shape B + gamma_b = solver_a.solve(body_b)["gamma_distribution"] + tables_b = solver_a.casadi_polar_table_inputs() + assert np.max(np.abs(tables_b["cl_table"] - tables_a["cl_table"])) > 1e-3 # the polars differ + + assert solver_a._casadi_newton_function(table_inputs=True) is fn + assert solver_a.casadi_panel_force_function(table_inputs=True) is force_fn + assert len(solver_a._casadi_table_input_cache) == cache_size + + fresh_b, _ = _solved(_v3("cosine")) + _assert_outputs_equal( + fn(**_newton_inputs(solver_a, gamma_b), **tables_b), + fresh_b._casadi_newton_function()(**_newton_inputs(solver_a, gamma_b)), + ) + _assert_outputs_equal( + force_fn(**_force_inputs(solver_a, gamma_b), **tables_b), + fresh_b.casadi_panel_force_function()(**_force_inputs(solver_a, gamma_b)), + ) + + +def test_omitted_tables_give_nan_not_zero(): + solver, gamma = _solved(_v3("uniform")) + out = solver._casadi_newton_function(table_inputs=True)(**_newton_inputs(solver, gamma)) + assert np.all(np.isnan(np.array(out["h"], dtype=float))) + out = solver.casadi_panel_force_function(table_inputs=True)(**_force_inputs(solver, gamma)) + assert np.all(np.isnan(np.array(out["force"], dtype=float))) + + +def test_table_input_functions_expand_to_sx(): + solver, gamma = _solved(_v3("uniform"), is_with_artificial_viscosity=True) + tables = solver.casadi_polar_table_inputs() + for fn, args in ( + (solver._casadi_newton_function(table_inputs=True), _newton_inputs(solver, gamma)), + (solver.casadi_panel_force_function(table_inputs=True), _force_inputs(solver, gamma)), + ): + _assert_outputs_equal(fn.expand()(**args, **tables), fn(**args, **tables)) + + +def test_newton_loop_is_unchanged_and_does_not_build_table_functions(): + solver, _ = _solved(_v3("uniform")) + assert solver._casadi_table_input_cache == {} + + +def test_unsupported_cases_raise_clearly(): + solver, _ = _solved(_v3("uniform"), polar_interpolation="bspline") + with pytest.raises(ValueError, match="polar_interpolation='linear'"): + solver._casadi_newton_function(table_inputs=True) + with pytest.raises(ValueError, match="polar_interpolation='linear'"): + solver.casadi_polar_table_inputs()