diff --git a/.gitignore b/.gitignore index d7c105ab..16f91689 100644 --- a/.gitignore +++ b/.gitignore @@ -59,3 +59,6 @@ lib/**/Manifest.toml # Generated test geometries (regenerated on demand by ram_air_matrix_wing) test/generated/ + +# Generated aero tables and run logs +*.arrow diff --git a/CHANGELOG.md b/CHANGELOG.md index 038f8be6..5a97726d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,23 +3,65 @@ ## VortexStepMethod v4.2.0 2026-08-25 ### Added + +- `SAMPLED` aero model: cl/cd/cm sampled at ascending angles of attack per + panel, interpolated between them and held flat past either end. This is what a + live polar source writes every solve; `set_sampled_polar!` rewrites a panel's + knots and values in place. It spans only the angles it was sampled over, so it + ignores `delta` and is skipped by the stall-angle scan — but inside that range + it represents a stall, which is what a local fit cannot do. +- Live in-memory polars in `AirfoilAero`: `KulfanBasis` and `deform_kulfan` + deform a fixed Kulfan fit analytically (a matvec against a constant CST basis, + never a refit, which is non-unique); `control_point_deflection` resamples a + deflection given at arbitrary chord fractions — beam nodes, membrane nodes — + onto that basis. `LivePolars` and `refresh_live_polars!` then evaluate + NeuralFoil on a grid of angles per panel in one batched forward pass and write + those values in as each panel's `SAMPLED` polar, so a chordwise deformation + reaches the aerodynamics as a shape change rather than a flap angle. The grid + moves with the panel — `LivePolarSettings` carries the offsets off its current + angle of attack — and `polar_drift` reports how far the solve has left it, + past which the polar holds its last sampled value instead of extrapolating. +- Live surface pressure in `AirfoilAero`: `refresh_live_pressure!` regenerates + every panel's `Cp` from its deformed shape in one batched forward pass at the + converged angle of attack, resampled onto the panel's own contour nodes by + `contour_pressure`; `live_surface_friction!` gives the matching skin friction + from the flat-plate closure at the panel's live Reynolds. A caller that + spreads panel forces over a surface can now make the pattern track the + deformation the forces already do. `contour_shape_matrix` and + `live_shape_offset!` move the contour itself by the same camber increment the + airfoil was deformed by, so the normals and segment areas a traction is built + from are the deformed section's rather than the reference one's — the chord + fractions never move, only the offsets, so a surface→point map built on that + contour stays valid. `decode_surface_pressure` and `contour_pressure` also + back the `NeuralFoilSolver` contour assembly, so the live path and the offline + table generator are one implementation. +- `prepare_inputs` accepts one Kulfan shape per case, so a whole wing goes + through NeuralFoil in a single forward pass. +- A panel carries the deformed airfoil its polar was generated from as + `live_shape`, and the Makie lofted-airfoil skin draws that object rather than + the tabulated `delta` contour when it is there. It is the shape the solver was + handed, not a re-derivation, so a deformation bug is visible in the picture. + `KulfanParameters` moved from `AirfoilAero` into the core module for it; it is + still exported from both. - Shared panel aerodynamics (`src/panel_aerodynamics.jl`): the per-panel physics - written once as pure, branch-free, number-type-generic functions of the section - geometry and the flow — `panel_axes`, `panel_inflow`, `panel_force_directions`, - `panel_loads` and the small helpers around them. `update_panel_properties!`, - `init_pos!`, `calc_forces!` and `calculate_results` now call them instead of - spelling the algebra out three times, and `SymbolicAWEModels` traces the same - functions with symbolic arguments to build its equations, so the two packages - can no longer drift apart. Panel geometry is unchanged bit for bit; forces, - moments and coefficients agree to within 1 ulp, the products having been - reassociated. `calc_forces!` stays zero-allocation. + written once as pure, branch-free, number-type-generic functions of the + section geometry and the flow — `panel_axes`, `panel_inflow`, + `panel_force_directions`, `panel_loads` and the small helpers around them. + `update_panel_properties!`, `init_pos!`, `calc_forces!` and + `calculate_results` now call them instead of spelling the algebra out three + times, and `SymbolicAWEModels` traces the same functions with symbolic + arguments to build its equations, so the two packages can no longer drift + apart. Panel geometry is unchanged bit for bit; forces, moments and + coefficients agree to within 1 ulp, the products having been reassociated. + `calc_forces!` stays zero-allocation. - `effective_alpha` and the `deficiency` argument of `panel_inflow` carry an unsteady lag (a Wagner indicial deficiency) into the angle the polars are read at, while the geometric angle still turns the force. Unused by the solver, - which has no unsteady state; it is shared so a symbolic consumer that does have - one reads the same definition. + which has no unsteady state; it is shared so a symbolic consumer that does + have one reads the same definition. ### Changed + - `section_pitch_rate` gained a three-argument form taking the trailing minus leading edge apparent wind directly. The four-argument form is unchanged. - Norms inside the shared panel aerodynamics are floored (`smooth_norm`) rather @@ -29,37 +71,43 @@ ## VortexStepMethod v4.1.2 2026-08-22 ### Fixed + - Billowing now takes its rotation axis from the section pair itself instead of projecting the leading-edge span vector on `spanwise_direction`. Sections run - `+y` to `-y`, so the pair already fixes the sign. The projection was decided by - 3 % of the vector on a C-shaped kite's tip panels — 11 mm of span between the - outermost sections — so a deforming tip could flip that one panel's billow. + `+y` to `-y`, so the pair already fixes the sign. The projection was decided + by 3 % of the vector on a C-shaped kite's tip panels — 11 mm of span between + the outermost sections — so a deforming tip could flip that one panel's + billow. ## VortexStepMethod v4.1.1 2026-08-19 ### Fixed + - Section twist now rotates the chord fully about the spanwise axis. The axial Rodrigues term was missing, so a swept or dihedral section's along-span chord - component was scaled by `cos(theta)`: the chord shortened and tilted out of the - twist plane, growing as `theta^2`. + component was scaled by `cos(theta)`: the chord shortened and tilted out of + the twist plane, growing as `theta^2`. ## VortexStepMethod v4.1.0 2026-08-17 ### Added -- `table_format` keyword on `write_section_aero`, `generate_airfoils`, `obj_to_yaml` - and `surfplan_to_aero_yaml`: `:csv` (default, readable) or `:arrow` (binary, ~40× - faster to load and 2.6× smaller). `read_section_aero` detects the format from the - file suffix, so a geometry YAML can reference either. + +- `table_format` keyword on `write_section_aero`, `generate_airfoils`, + `obj_to_yaml` and `surfplan_to_aero_yaml`: `:csv` (default, readable) or + `:arrow` (binary, ~40× faster to load and 2.6× smaller). `read_section_aero` + detects the format from the file suffix, so a geometry YAML can reference + either. - `geometry_path` keyword on `obj_to_yaml`, naming the geometry YAML itself instead of always writing `output_dir/geometry.yaml`. Point it outside the - table directory and the emitted table references carry the path from the YAML's - directory to `output_dir`, which is what the geometry loader resolves them - against — so a generated dataset can keep its bulk in a subdirectory while the - geometry sits with the hand-written ones. -- `convert_node_table` and `write_node_rows` rewrite a per-node table in the format - the destination suffix names. `obj_to_yaml` migrates an existing dataset with - them when `table_format` differs from what the directory holds, so a dataset - changes format without re-running the airfoil solver that produced it. + table directory and the emitted table references carry the path from the + YAML's directory to `output_dir`, which is what the geometry loader resolves + them against — so a generated dataset can keep its bulk in a subdirectory + while the geometry sits with the hand-written ones. +- `convert_node_table` and `write_node_rows` rewrite a per-node table in the + format the destination suffix names. `obj_to_yaml` migrates an existing + dataset with them when `table_format` differs from what the directory holds, + so a dataset changes format without re-running the airfoil solver that + produced it. - `flow_curvature` solver setting (default `false`). When enabled, each section gets the thin-airfoil pitch-rate moment increment `Δcm = -(π/4) q̂` with `q̂ = q c / (2 v_rel)`. A section rotating about its own spanwise axis sees an @@ -68,131 +116,146 @@ point cannot represent. The lift response to `q` was already exact because the inflow is sampled at the three-quarter-chord point, so only the moment was missing. -- `pitch_rate_dist` field on `BodyAerodynamics`, holding the `q` above per panel. - `set_va!(body_aero, va, omega)` fills it by projecting `omega` onto every - panel's `y_airf`, so panels at different dihedral see different rates from one - body rate. The distributed `set_va!(body_aero, va_distribution; - pitch_rate_dist)` takes it directly, so twist and flapping rates of a deforming - wing — which no single body rate can express — reach the moment. Omitting the - keyword zeroes it rather than reusing a stale `omega`, which the distributed - form never sets. -- `section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord)` builds - one entry of that distribution from a section's edge velocities, and reduces to - `ω ⋅ y_airf` for rigid motion. -- `bin/update_default_manifests` regenerates both checked-in default manifests in - one command, and `bin/install` takes `+X.Y` / `--version X.Y` to pick a Julia - channel without prompting. Selecting a version no longer moves the juliaup - default, so regenerating the 1.11 manifest leaves the shell on whatever channel - it was using. +- `pitch_rate_dist` field on `BodyAerodynamics`, holding the `q` above per + panel. `set_va!(body_aero, va, omega)` fills it by projecting `omega` onto + every panel's `y_airf`, so panels at different dihedral see different rates + from one body rate. The distributed + `set_va!(body_aero, va_distribution; + pitch_rate_dist)` takes it directly, so + twist and flapping rates of a deforming wing — which no single body rate can + express — reach the moment. Omitting the keyword zeroes it rather than reusing + a stale `omega`, which the distributed form never sets. +- `section_pitch_rate(velocity_leading, velocity_trailing, z_airf, chord)` + builds one entry of that distribution from a section's edge velocities, and + reduces to `ω ⋅ y_airf` for rigid motion. +- `bin/update_default_manifests` regenerates both checked-in default manifests + in one command, and `bin/install` takes `+X.Y` / `--version X.Y` to pick a + Julia channel without prompting. Selecting a version no longer moves the + juliaup default, so regenerating the 1.11 manifest leaves the shell on + whatever channel it was using. ### Fixed + - The `Solver` docstring quoted the `SolverSettings` defaults for - `type_initial_gamma_distribution` and `core_radius_fraction` (`ELLIPTIC`, `1e-20`) - instead of its own (`ZEROS`, `0.05`), and never said what `core_radius_fraction` - measures. It now documents the `Solver` defaults and cites Damiani et al. (2019) for - the 0.05 cut-off. + `type_initial_gamma_distribution` and `core_radius_fraction` (`ELLIPTIC`, + `1e-20`) instead of its own (`ZEROS`, `0.05`), and never said what + `core_radius_fraction` measures. It now documents the `Solver` defaults and + cites Damiani et al. (2019) for the 0.05 cut-off. - A remesh under `use_prior_polar` no longer resamples the refined sections' `SectionAero` surface tables down to whatever unrefined sections survive it. - `compute_refined_section_interpolation!` reblended contour, `cp` and `cf` from the - unrefined sections unconditionally while `aero_data` was preserved, so a wing rebuilt - onto fewer structural stations kept full-resolution polars but lost the surface tables - pressure integration reads. The reblend is now skipped when the polars are preserved - and the refined sections already carry tables. + `compute_refined_section_interpolation!` reblended contour, `cp` and `cf` from + the unrefined sections unconditionally while `aero_data` was preserved, so a + wing rebuilt onto fewer structural stations kept full-resolution polars but + lost the surface tables pressure integration reads. The reblend is now skipped + when the polars are preserved and the refined sections already carry tables. - The angle-of-attack correction (`correct_aoa=true`) no longer overwrites - `body_aero.AIC` with the aerodynamic-centre (LLT) matrix. It builds that matrix in - its own `AIC_aero_center` field, so `AIC` keeps holding the control-point matrix the - circulation was solved against, which is what anything reading it after a solve - expects. Costs a second `n_panels × n_panels × 3` buffer. -- `obj_to_yaml` no longer places sections on a wingtip that has closed to a point. - `station_indices` spreads its targets over the stations that still have a chord, - so the outermost section lands on the last sliceable one. A V3 mesh sliced with - the default `wingtip_distance` used to put a zero-chord section at each tip, - whose polar was `NaN` and took the whole solve with it; working around it meant - guessing a `wingtip_distance` large enough to skip past the tip. That workaround - is no longer the default: `wingtip_distance` is now `0.0`, an inset on top of the - trim for meshes whose slices just short of the tip are still too thin to analyse. -- Wing sections are normalized to `+y` to `-y` order on load (`normalize_span_order!`), - and by `refine!` for wings built through `add_section!`. Panel `y_airf` and `z_airf` - follow the order sections are stored in, so a geometry file written the other way - round inverted every panel normal, and a wing whose sections were replaced after its - panels were built (a structural remesh) inverted them mid-run, one panel at a time as - each crossed `spanwise_direction`. `obj_to_yaml` and `surfplan_to_aero_yaml` emit that + `body_aero.AIC` with the aerodynamic-centre (LLT) matrix. It builds that + matrix in its own `AIC_aero_center` field, so `AIC` keeps holding the + control-point matrix the circulation was solved against, which is what + anything reading it after a solve expects. Costs a second + `n_panels × n_panels × 3` buffer. +- `obj_to_yaml` no longer places sections on a wingtip that has closed to a + point. `station_indices` spreads its targets over the stations that still have + a chord, so the outermost section lands on the last sliceable one. A V3 mesh + sliced with the default `wingtip_distance` used to put a zero-chord section at + each tip, whose polar was `NaN` and took the whole solve with it; working + around it meant guessing a `wingtip_distance` large enough to skip past the + tip. That workaround is no longer the default: `wingtip_distance` is now + `0.0`, an inset on top of the trim for meshes whose slices just short of the + tip are still too thin to analyse. +- Wing sections are normalized to `+y` to `-y` order on load + (`normalize_span_order!`), and by `refine!` for wings built through + `add_section!`. Panel `y_airf` and `z_airf` follow the order sections are + stored in, so a geometry file written the other way round inverted every panel + normal, and a wing whose sections were replaced after its panels were built (a + structural remesh) inverted them mid-run, one panel at a time as each crossed + `spanwise_direction`. `obj_to_yaml` and `surfplan_to_aero_yaml` emit that order too; files of either order keep loading the same. -- Spanwise distribution plots put `+y` on the left, matching that order and the kite - seen from the front. -- The lofted airfoil skin in `plot_geometry` draws each section's contour on its own - panel edge (#256). The edge was picked from the wing's span order rather than from - the panel index, so a wing whose sections ran `-y` to `+y` got every rib drawn at its - neighbour's station, one tip bare and the other doubled, while the panels themselves - rendered correctly. +- Spanwise distribution plots put `+y` on the left, matching that order and the + kite seen from the front. +- The lofted airfoil skin in `plot_geometry` draws each section's contour on its + own panel edge (#256). The edge was picked from the wing's span order rather + than from the panel index, so a wing whose sections ran `-y` to `+y` got every + rib drawn at its neighbour's station, one tip bare and the other doubled, + while the panels themselves rendered correctly. ### Changed -- `SolverSettings` now defaults to the same values as `Solver`: `core_radius_fraction` - `1e-20` → `0.05` and `type_initial_gamma_distribution` `ELLIPTIC` → `ZEROS`. The 0.05 - bound vortex core cut-off follows Damiani et al. (2019), "A Vortex Step Method for - Nonlinear Airfoil Polar Data as Implemented in KiteAeroDyn", and matches the upstream - `awegroup/Vortex-Step-Method` default; at `1e-20` the Biot-Savart singularity guard - never engaged. Coefficients are unchanged for well-separated geometry, since the guard - only triggers where a control point falls within 5% of a filament length of a bound - vortex, and every settings file shipped in `data/` sets both keys explicitly. + +- `SolverSettings` now defaults to the same values as `Solver`: + `core_radius_fraction` `1e-20` → `0.05` and `type_initial_gamma_distribution` + `ELLIPTIC` → `ZEROS`. The 0.05 bound vortex core cut-off follows Damiani et + al. (2019), "A Vortex Step Method for Nonlinear Airfoil Polar Data as + Implemented in KiteAeroDyn", and matches the upstream + `awegroup/Vortex-Step-Method` default; at `1e-20` the Biot-Savart singularity + guard never engaged. Coefficients are unchanged for well-separated geometry, + since the guard only triggers where a control point falls within 5% of a + filament length of a bound vortex, and every settings file shipped in `data/` + sets both keys explicitly. - `BodyAerodynamics.AIC` is stored as `(n_panels, n_panels, 3)` instead of - `(3, n_panels, n_panels)`, so each component slice `AIC[:, :, k]` is contiguous and - the induced-velocity products reach BLAS `gemv` instead of the generic fallback. - `solve!` is 3.1–5.3× faster (n=120, VSM: 26.1 ms → 4.9 ms inviscid, 21.1 ms → - 6.0 ms with polars). Code reading `AIC[k, i, j]` must become `AIC[i, j, k]`. -- The filament induced-velocity kernels reuse the `r0`/`length` each `BoundFilament` - already stores, take the core-radius cutoff from `|r1.r0|/|r0|` without forming the - perpendicular vector, and defer the cross products that only one branch reads. Output - is bit-identical; `solve!` is a further 1.47-1.55x faster. -- `read_node_table` parses into a preallocated matrix instead of `reduce(vcat, …)` - over a generator, which was quadratic in the row count: ~21× faster on a 16 MB - surface table (2.49 s → 0.12 s), benefiting every existing dataset. + `(3, n_panels, n_panels)`, so each component slice `AIC[:, :, k]` is + contiguous and the induced-velocity products reach BLAS `gemv` instead of the + generic fallback. `solve!` is 3.1–5.3× faster (n=120, VSM: 26.1 ms → 4.9 ms + inviscid, 21.1 ms → 6.0 ms with polars). Code reading `AIC[k, i, j]` must + become `AIC[i, j, k]`. +- The filament induced-velocity kernels reuse the `r0`/`length` each + `BoundFilament` already stores, take the core-radius cutoff from + `|r1.r0|/|r0|` without forming the perpendicular vector, and defer the cross + products that only one branch reads. Output is bit-identical; `solve!` is a + further 1.47-1.55x faster. +- `read_node_table` parses into a preallocated matrix instead of + `reduce(vcat, …)` over a generator, which was quadratic in the row count: ~21× + faster on a 16 MB surface table (2.49 s → 0.12 s), benefiting every existing + dataset. - `is_show=true` draws into a window named after the plot title instead of into whichever window the backend last used, so a script showing several plots gets one window each and re-running it redraws them in place. `show_plot` takes the window `name` as a keyword. -- `examples/V3_kite.jl` builds its wing from `VSMSettings` and adds a second sweep on - polars generated from `V3_25.obj` with NeuralFoil (`NEURALFOIL`, on by default), so - the plots compare CAD-derived polars against the checked-in CFD tables. -- `bin/run_julia` starts Julia with `JULIA_NUM_THREADS=auto` unless the environment - already sets it. +- `examples/V3_kite.jl` builds its wing from `VSMSettings` and adds a second + sweep on polars generated from `V3_25.obj` with NeuralFoil (`NEURALFOIL`, on + by default), so the plots compare CAD-derived polars against the checked-in + CFD tables. +- `bin/run_julia` starts Julia with `JULIA_NUM_THREADS=auto` unless the + environment already sets it. ## VortexStepMethod v4.0.0 2026-08-03 ### Added -- `NeuralFoil`-based airfoil polar generation via the new `AirfoilAero` submodule - (`NeuralFoilSolver`, `XFoilSolver`, `analyze_section`, `analyze_sweep`, - `fit_kulfan_parameters`, `shrink_wrap`, `ShrinkWrap`) + +- `NeuralFoil`-based airfoil polar generation via the new `AirfoilAero` + submodule (`NeuralFoilSolver`, `XFoilSolver`, `analyze_section`, + `analyze_sweep`, `fit_kulfan_parameters`, `shrink_wrap`, `ShrinkWrap`) - `ObjAdapter` submodule: converts a 3D wing `.obj` mesh to the native YAML/CSV - geometry format (`obj_to_yaml`, `perpendicular_sections`, - `write_yaml`, `plot_slices_3d`, `plot_airfoils`) -- `SurfplanAdapter` submodule: `surfplan_to_aero_yaml` turns a SurfplanAdapter aero - export into the native pressure-ready YAML/CSV geometry (shared `generate_airfoils` - core with `ObjAdapter`) -- `ObjWing(obj_path[, dat_path]; Re, n_panels, aero_solver, remake, ...)` convenience - constructor restored for backward compatibility — internally calls + geometry format (`obj_to_yaml`, `perpendicular_sections`, `write_yaml`, + `plot_slices_3d`, `plot_airfoils`) +- `SurfplanAdapter` submodule: `surfplan_to_aero_yaml` turns a SurfplanAdapter + aero export into the native pressure-ready YAML/CSV geometry (shared + `generate_airfoils` core with `ObjAdapter`) +- `ObjWing(obj_path[, dat_path]; Re, n_panels, aero_solver, remake, ...)` + convenience constructor restored for backward compatibility — internally calls `ObjAdapter.obj_to_yaml` then `Wing`; `aero_solver` selects the polar backend (default `NeuralFoilSolver()`; pass `XFoilSolver()` for old behavior); - `remake=false` (default) reuses an existing `geometry.yaml` in `output_dir` to skip - expensive polar generation when only `n_panels` changes; set `remake=true` to force - regeneration -- Per-section surface aero tables: `SectionAero` (full contour + surface pressure - `cp` + skin friction `cf` per node over an `(α, δ)` grid), `section_surface`, - `read_section_aero` / `write_section_aero`; surface aero propagates through - `Section` and `Wing` and is spanwise-interpolated during `refine!` + `remake=false` (default) reuses an existing `geometry.yaml` in `output_dir` to + skip expensive polar generation when only `n_panels` changes; set + `remake=true` to force regeneration +- Per-section surface aero tables: `SectionAero` (full contour + surface + pressure `cp` + skin friction `cf` per node over an `(α, δ)` grid), + `section_surface`, `read_section_aero` / `write_section_aero`; surface aero + propagates through `Section` and `Wing` and is spanwise-interpolated during + `refine!` - `POLY` aero model for polynomial cl/cd/cm (exported) -- `lei_poly_coeffs(tube_diameter, camber)` (exported from `AirfoilAero`) returning - the Breukels α-polynomial cl/cd/cm coefficients +- `lei_poly_coeffs(tube_diameter, camber)` (exported from `AirfoilAero`) + returning the Breukels α-polynomial cl/cd/cm coefficients - Li/Gaunaa spanwise artificial viscosity (Li, Gaunaa, Pirrung & Lønbæk, TORQUE 2026) for the LOOP solver, stabilizing post-stall circulation distributions that otherwise develop non-physical sawtooth oscillations; opt-in via `is_with_artificial_viscosity` (default `false`) and `artificial_viscosity_factor` (default `0.035`) on the solver settings -- `crease_frac` wing setting (chordwise flap-hinge fraction, default `0.75`) - for drawing the δ-deflected plate/skin, readable from YAML +- `crease_frac` wing setting (chordwise flap-hinge fraction, default `0.75`) for + drawing the δ-deflected plate/skin, readable from YAML - `examples/V3_neuralfoil.jl` and `examples/obj_to_yaml_kite.jl` ### Changed + - BREAKING - `ObjWing` polar generation now uses NeuralFoil (via `ObjAdapter.obj_to_yaml`) instead of XFoil + a user-supplied `.dat` file. The constructor signature is backward-compatible (`dat_path` is accepted and @@ -201,26 +264,29 @@ `aero_solver=XFoilSolver()` to `ObjAdapter.obj_to_yaml` and call `Wing` directly. - OBJ-based wings are now built via `ObjAdapter.obj_to_yaml` + `Wing(yaml_path)` - instead of the old `ObjWing` pipeline (XFoil + single `.dat` file); - `ObjWing` is kept as a shim that accepts but ignores `dat_path` + instead of the old `ObjWing` pipeline (XFoil + single `.dat` file); `ObjWing` + is kept as a shim that accepts but ignores `dat_path` - BREAKING - `LEI_AIRFOIL_BREUKELS` is now a deprecated alias of `POLY`, and the built-in Breukels regression no longer runs at solve time. Sections must carry the α-polynomial coefficients directly instead of `(tube_diameter, camber)`; compute them up front with `lei_poly_coeffs(tube_diameter, camber)` - `Section` and `add_section!` accept an optional `section_aero` argument -- BREAKING - plotting is now Makie-only; the `VortexStepMethodMakieExt` extension - loads once a Makie backend and - [`MakieControlPlots`](https://github.com/OpenSourceAWE/MakieControlPlots.jl) are - available, so plotting code must now load `MakieControlPlots` (and drop - `set_plot_backend!`); `plot_section_polars` is rendered through `MakieControlPlots` +- BREAKING - plotting is now Makie-only; the `VortexStepMethodMakieExt` + extension loads once a Makie backend and + [`MakieControlPlots`](https://github.com/OpenSourceAWE/MakieControlPlots.jl) + are available, so plotting code must now load `MakieControlPlots` (and drop + `set_plot_backend!`); `plot_section_polars` is rendered through + `MakieControlPlots` ### Removed + - `ObjWing` as a standalone pipeline (replaced by `ObjAdapter`); the name is re-exported as a compatibility wrapper - `auto_rotation` helper (internal, removed from public API) -- `PanelGroupingMethod` enum (already removed in v3.0.0 — stale docs entry cleaned up) -- ControlPlots test run removed from CI (`plot-controlplots` arg) due to - a `libraqm`/HarfBuzz symbol conflict in the GitHub Actions environment +- `PanelGroupingMethod` enum (already removed in v3.0.0 — stale docs entry + cleaned up) +- ControlPlots test run removed from CI (`plot-controlplots` arg) due to a + `libraqm`/HarfBuzz symbol conflict in the GitHub Actions environment - BREAKING - the `ControlPlots` plotting extension, `examples_cp/`, and the `PythonCall`/Matplotlib setup (`bin/install_controlplots`, CondaPkg `LocalPreferences` defaults) @@ -232,6 +298,7 @@ ## VortexStepMethod v3.3.6 2026-06-13 ### Added + - `calc_forces!` and `solve_base!` (both exported): `solve!` is now `solve_base!` followed by `calc_forces!`, so a frozen circulation can be mapped to forces without re-running the nonlinear gamma solve (#245) @@ -240,10 +307,12 @@ kept as a thin wrapper (#246) ### Changed -- `calc_forces!` is now allocation-free in the per-step hot path - (preallocated `panel_area_dist` and `unrefined_count_dist` buffers) (#245) + +- `calc_forces!` is now allocation-free in the per-step hot path (preallocated + `panel_area_dist` and `unrefined_count_dist` buffers) (#245) ### Fixed + - 3D polar plotting (#245) - flaky Aqua `persistent_tasks` test now actually disabled via `persistent_tasks=false` (`()` did not disable it) (#246) @@ -251,19 +320,22 @@ ## VortexStepMethod v3.3.5 2026-06-05 ### Added + - `moment_coeff_unrefined_dist` field in `VSMSolution`: the summed `moment_frac`-referenced pitching-moment coefficient per unrefined section [-] ## VortexStepMethod v3.3.4 2026-05-31 ### Added -- `PlotBackend`, `MakieBackend`, `ControlPlotsBackend`, and - `set_plot_backend!` so applications can explicitly choose which plotting - extension the backend-agnostic plotting API should use -- `PythonCall` added as a weak dependency to support the `ControlPlots` - backend with PythonPlot + +- `PlotBackend`, `MakieBackend`, `ControlPlotsBackend`, and `set_plot_backend!` + so applications can explicitly choose which plotting extension the + backend-agnostic plotting API should use +- `PythonCall` added as a weak dependency to support the `ControlPlots` backend + with PythonPlot ### Changed + - backend-agnostic plotting wrappers now route through the active plotting backend, and each plotting extension initializes itself as the default only when no backend has been selected yet @@ -271,6 +343,7 @@ - improved `bin/install` and `bin/install_controlplots` scripts ### Fixed + - corrected projection onto core radius in `velocity_3D_bound_vortex!` and semi-infinite trailing vortex projection so that the radial direction is always measured from the filament axis, not from the origin (#241) @@ -281,6 +354,7 @@ ## VortexStepMethod v3.3.3 2026-05-21 ### Fixed + - `MakieExt` and `ControlPlotsExt` no longer both define `VortexStepMethod.plot_geometry` for the same type, resolving a method ambiguity when both extensions were loaded (#236) @@ -290,33 +364,39 @@ ## VortexStepMethod v3.3.2 2026-05-18 ### Changed + - use 2-arg version of atan to avoid possible NaN ## VortexStepMethod v3.3.1 2026-05-13 ### Changed + - `unrefined_deform!` linearly interpolates twist and TE deflection between unrefined sections and rotates each refined section about the average of its adjacent local airfoil normals (#234) ### Fixed -- `smooth_sqrt` in the solver hot loop keeps gradients defined at zero - velocity magnitude + +- `smooth_sqrt` in the solver hot loop keeps gradients defined at zero velocity + magnitude ## VortexStepMethod v3.3.0 2026-05-05 ### Added + - `ForwardDiff` compatibility, used by default in `linearize` (#232) - `backend` keyword argument for `linearize` - example `linearize_check.jl` comparing FiniteDiff and ForwardDiff tangents ### Changed + - core structs are parameterized on the scalar type `T` so dual numbers can propagate through them; public constructors are unchanged ## VortexStepMethod v3.2.0 2026-05-02 ### Added + - support for both CairoMakie and GLMakie in the example `menu()` via new `CairoMakie.activate()` / `GLMakie.activate()` entries (#222) - plots are saved as PDF when CairoMakie is active @@ -324,18 +404,21 @@ generated plots ### Changed + - examples write plots into the shared `output` folder instead of the working directory - example plot file names are sanitized: spaces replaced with `_` and `%` replaced with `pct` ### Fixed + - NONLIN solver no longer returns stale `gamma` on repeated `solve!` calls (#228) ## VortexStepMethod v3.1.3 2026-04-23 ### Fixed + - bug in `linearize` where `body_aero.va` was used instead of `body_aero._va`, causing incorrect initial velocity storage (#227) - `set_va!` no longer overwrites `_va` with a computed reference velocity when @@ -344,46 +427,63 @@ - linearize now correctly preserves and restores `omega` across perturbations ## VortexStepMethod v3.1.2 2025-04-20 + ### Added + - add back compat entry v2 for SciMLBase ## VortexStepMethod v3.1.1 2025-04-20 + ### Added + - add back compat entry v3 for RecursiveArrayTools ## VortexStepMethod v3.1.0 2025-04-19 ### Breaking + - `billowing_angle` replaced by `billowing_percentage` on `Wing` and `WingSettings` (percentage of arc length, not radians) - `billowing_angle_from_percentage()` removed - `BILLOWING` distribution now uses sinusoidal rotation instead of circular arc ### Added + - `billowing.jl` example comparing flat vs billowed V3 kite - Coarse V3 kite geometry, settings, and combined CFD polar data - `cl_over_cd` keyword for `plot_polars` and `plot_combined_analysis` -- the function `menu_cp()` can now be used to run the examples with the ControlPlots backend -- the script `bin/install` can and should be used to instantiate the project and all sub-projects after git checkout -- the files `Manifest-v1.11.toml.default` and `Manifest-v12.toml.default` for enhanced reproducibility -- the scripts `bin/jetls` and `bin/jetls_examples` for running a static check of the source code -- the script `test_bench.jl` for benchmarking the refinement method and measuring allocations +- the function `menu_cp()` can now be used to run the examples with the + ControlPlots backend +- the script `bin/install` can and should be used to instantiate the project and + all sub-projects after git checkout +- the files `Manifest-v1.11.toml.default` and `Manifest-v12.toml.default` for + enhanced reproducibility +- the scripts `bin/jetls` and `bin/jetls_examples` for running a static check of + the source code +- the script `test_bench.jl` for benchmarking the refinement method and + measuring allocations ### Changed + - Plot legends moved to shared horizontal legend at bottom of grid layouts - the script bin/run_julia now can be called with a script name as parameter -- fixed all JETLS warnings in the source code for improved performance and stability +- fixed all JETLS warnings in the source code for improved performance and + stability ### Fixed errors -- Domain error in elliptical gamma distribution when control points lie - outside the nominal span envelope + +- Domain error in elliptical gamma distribution when control points lie outside + the nominal span envelope ## VortexStepMethod v3.0.1 2025-04-04 ### Changed + - the file `CITATION.cff` - compat entry for RecursiveArrayTools + ### Added + - the file `.zenodo.json` ## VortexStepMethod v3.0.0 diff --git a/docs/src/functions.md b/docs/src/functions.md index 6562b163..dc7fdcf2 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -41,6 +41,19 @@ deform_section analyze_section analyze_sweep neuralfoil_aero +deform_kulfan +chord_residual +chord_line +control_point_deflection +panel_kulfan_parameters +refresh_live_polars! +deform_live_shapes! +apply_live_shapes! +polar_drift +refresh_live_pressure! +live_surface_friction! +contour_shape_matrix +live_shape_offset! generate_aero_matrices generate_polar_from_coordinates generate_polar_from_dat diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index b48ed3e8..e657184c 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -98,6 +98,10 @@ update_non_deformed_sections! ### Aerodynamic data and Cp ```@docs calculate_new_aero_data +set_sampled_polar! +decode_surface_pressure +contour_pressure +sampled_value assemble_polar_matrix load_matrix_polar_data read_aero_matrix @@ -153,8 +157,11 @@ smooth_turning! load_neuralfoil_model neuralfoil_section neuralfoil_fused_output +fused_output +decode_coefficients nn_forward prepare_inputs +fill_case_input! flip_inputs flip_outputs squared_mahalanobis_distance @@ -221,6 +228,7 @@ map_airfoil_3d fitted_airfoil_3d generated_slices airfoil_skin_geometry +panel_contour panel_normal plate_hinge_local panel_plate_geometry diff --git a/docs/src/private_types.md b/docs/src/private_types.md index 5fca5ebe..ff63c1fc 100644 --- a/docs/src/private_types.md +++ b/docs/src/private_types.md @@ -7,6 +7,7 @@ CurrentModule = VortexStepMethod ### Wing Geometry, Panel and Aerodynamics ```@docs Panel +KulfanParameters PanelProperties Filament BoundFilament @@ -23,7 +24,9 @@ LEI_AIRFOIL_BREUKELS CurrentModule = VortexStepMethod.AirfoilAero ``` ```@docs -KulfanParameters +KulfanBasis +LivePolarSettings +LivePolars NeuralFoilModel NeuralFoilResult ``` diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index 1356df73..ddd010e2 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/ext/VortexStepMethodMakieExt.jl @@ -146,15 +146,34 @@ function Makie.plot!(ax, panel::VortexStepMethod.Panel; color=(:red, 0.2), R_b_w return plots end +""" + panel_contour(panel, section) -> (x, y) + +The airfoil contour to draw a panel with: coordinates of its `live_shape` when a live +polar source has put one there, else the section's tabulated contour at the panel's +`delta`. The live branch reads the stored shape object itself, so what is drawn is what +was flown. +""" +function panel_contour(panel, section) + shape = panel.live_shape + isnothing(shape) || return VortexStepMethod.AirfoilAero.kulfan_to_coordinates(shape) + return VortexStepMethod.section_surface(section.section_aero, 0.0, panel.delta)[1:2] +end + """ airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) -> (vertices, faces, ribs) -Lofted airfoil skin of a `BodyAerodynamics`: each section's deflected contour -(`section_surface` at the panel's `delta`) is fitted between the panel's `corner_points` -by a 2D similarity, TE pinned so a deflection bulges the fore body up. The skin reflects -`delta` only when the geometry carries per-`delta` slices (`obj_to_yaml` with a -`delta_range`); with δ=0-only data it renders undeflected, a deliberate cue that the -deflected slices are missing. Transformed to world by `R_b_w`/`T_b_w`. +Lofted airfoil skin of a `BodyAerodynamics`: each section's contour is fitted between +the panel's `corner_points` by a 2D similarity, TE pinned so a deflection bulges the +fore body up. + +The contour is the panel's own `live_shape` when it has one — the very +[`KulfanParameters`](@ref) object the live polar source deformed and handed to the +airfoil solver, not a re-derivation of it, so a deformation bug shows up in the picture +instead of being papered over. Otherwise it is `section_surface` at the panel's `delta`, +which reflects `delta` only when the geometry carries per-`delta` slices (`obj_to_yaml` +with a `delta_range`); with δ=0-only data it renders undeflected, a deliberate cue that +the deflected slices are missing. Transformed to world by `R_b_w`/`T_b_w`. `vertices`/`faces` triangulate the skin between consecutive equal-node sections; `ribs` is one closed contour polyline per section. Sections without contour data are skipped. """ @@ -169,10 +188,10 @@ function airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) n_panels = n - 1 n_panels < 1 && continue for (i, section) in enumerate(sections) - isnothing(section.section_aero) && continue panel_idx = panel_offset + min(i, n_panels) panel_idx <= length(body.panels) || continue panel = body.panels[panel_idx] + (isnothing(section.section_aero) && isnothing(panel.live_shape)) && continue corners = panel.corner_points corner1 = Point3f(corners[:, 1]); corner3 = Point3f(corners[:, 3]) corner2 = Point3f(corners[:, 2]); corner4 = Point3f(corners[:, 4]) @@ -183,8 +202,7 @@ function airfoil_skin_geometry(body; R_b_w=nothing, T_b_w=nothing) chord_len = norm(chord) chord_len < 1e-9 && continue up = panel_normal(panel) - xs, ys, _, _ = VortexStepMethod.section_surface(section.section_aero, - 0.0, panel.delta) + xs, ys = panel_contour(panel, section) le_i = argmin(xs) le_x = xs[le_i]; le_y = ys[le_i] te_x = 0.5 * (xs[1] + xs[end]); te_y = 0.5 * (ys[1] + ys[end]) diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 087a346a..62b7a859 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -36,7 +36,8 @@ export calculate_projected_area, calculate_span export MVec3 export LLT, Model, VSM -export AeroModel, INVISCID, POLY, LEI_AIRFOIL_BREUKELS, POLAR_MATRICES, POLAR_VECTORS +export AeroModel, INVISCID, POLY, LEI_AIRFOIL_BREUKELS, POLAR_MATRICES, POLAR_VECTORS, SAMPLED +export KulfanParameters export BILLOWING, COSINE, LINEAR, PanelDistribution, SPLIT_PROVIDED, UNCHANGED export ELLIPTIC, InitialGammaDistribution, ZEROS export FAILURE, FEASIBLE, INFEASIBLE, SolverStatus @@ -237,7 +238,7 @@ Enumeration of the implemented wing types. @enum WingType RECTANGULAR CURVED ELLIPTICAL """ - AeroModel `POLY` `POLAR_VECTORS` `POLAR_MATRICES` `INVISCID` + AeroModel `POLY` `POLAR_VECTORS` `POLAR_MATRICES` `INVISCID` `SAMPLED` Enumeration of the implemented aerodynamic models. See also: [AeroData](@ref) @@ -247,6 +248,11 @@ Enumeration of the implemented aerodynamic models. See also: [AeroData](@ref) - `POLAR_VECTORS`: Polar vectors as function of alpha (lookup tables with interpolation) - `POLAR_MATRICES`: Polar matrices as function of alpha and delta (lookup tables with interpolation) - INVISCID +- `SAMPLED`: cl/cd/cm sampled at ascending angles of attack per panel, interpolated + between them and held flat past either end. This is what a live polar source writes + each solve, see [`refresh_live_polars!`](@ref + VortexStepMethod.AirfoilAero.refresh_live_polars!). Samples are the polar, so a stall + knee inside the sampled range is represented rather than smoothed. Ignores `delta`. `LEI_AIRFOIL_BREUKELS` is a deprecated alias of `POLY`. @@ -257,6 +263,7 @@ where `alpha` is the angle of attack, `delta` is trailing edge angle. POLAR_VECTORS POLAR_MATRICES INVISCID + SAMPLED end """ @@ -339,6 +346,9 @@ Union of different definitions of the aerodynamic properties of a wing section. - (`alpha_range`, `cl_vector`, `cd_vector`, `cm_vector`) for `POLAR_VECTORS` - (`alpha_range`, `delta_range`, `cl_matrix`, `cd_matrix`, `cm_matrix`) for `POLAR_MATRICES` +`SAMPLED` carries no section-level data: it is written onto a panel at run time by a +live polar source, never read from a section. + where `alpha` is the angle of attack [rad], `delta` is trailing edge angle [rad], `cl` the lift coefficient, `cd` the drag coefficient and `cm` the pitching moment coefficient. The camber of a kite refers to the curvature of its airfoil shape. The camber is typically measured as the maximum distance diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index 531a5ad2..d02b1acd 100644 --- a/src/airfoil_aero/AirfoilAero.jl +++ b/src/airfoil_aero/AirfoilAero.jl @@ -8,11 +8,14 @@ using NPZ using Xfoil using Printf: @sprintf using ..VortexStepMethod: SectionAero, interpolate_matrix_nans!, delta_suffix, - write_node_rows + write_node_rows, section_surface, set_sampled_polar!, + KulfanParameters include("kulfan.jl") +include("deform.jl") include("shrink_wrap.jl") include("neuralfoil.jl") +include("live_polar.jl") include("poly.jl") include("airfoil_solvers/common.jl") include("airfoil_solvers/xfoil_solver.jl") @@ -29,6 +32,12 @@ export ShrinkWrap, shrink_wrap export fit_kulfan_parameters, kulfan_to_coordinates export NeuralFoilModel, NeuralFoilResult, load_neuralfoil_model export neuralfoil_aero, neuralfoil_section +export KulfanBasis, deform_kulfan, control_point_deflection +export chord_residual, chord_line +export LivePolarSettings, LivePolars, panel_kulfan_parameters +export refresh_live_polars!, refresh_live_pressure!, live_surface_friction! +export contour_shape_matrix, live_shape_offset! +export polar_drift export AbstractAirfoilSolver, XFoilSolver, NeuralFoilSolver export SectionSolution, DeformedSection, deform_section, analyze_section, analyze_sweep export create_2d_polars, generate_aero_matrices, generate_section_aero, lei_poly_coeffs diff --git a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl index 9a625600..6a807d39 100644 --- a/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl +++ b/src/airfoil_aero/airfoil_solvers/neuralfoil_solver.jl @@ -47,9 +47,8 @@ Assemble a full-contour [`SectionSolution`](@ref) for case `i` of a NeuralFoil s `(x, y)` split at the leading edge `le`, carrying the precomputed `cf`. """ function neuralfoil_contour_solution(alpha, res, i, x, y, le, cf) - up = linear_interpolation(res.x, res.cp_upper[:, i]; extrapolation_bc=Line()) - lo = linear_interpolation(res.x, res.cp_lower[:, i]; extrapolation_bc=Line()) - cp = [(k <= le ? up : lo)(clamp(x[k], 0.0, 1.0)) for k in eachindex(x)] + cp = contour_pressure(res.x, view(res.cp_upper, :, i), + view(res.cp_lower, :, i), x, le) return SectionSolution(alpha, res.cl[i], res.cd[i], res.cm[i], res.confidence[i], x, y, cp, cf) end diff --git a/src/airfoil_aero/deform.jl b/src/airfoil_aero/deform.jl new file mode 100644 index 00000000..c0f551fa --- /dev/null +++ b/src/airfoil_aero/deform.jl @@ -0,0 +1,145 @@ +""" + KulfanBasis(; n_stations=60, n_weights=8, ridge=1e-3) + +The fixed CST basis a shape deformation is projected onto: chord stations `x` in +`[0, 1]` and the ridge-regularised inverse `projection` of the class-times-Bernstein +matrix `C(x)·B(x)` those stations span. + +`ridge` is the Tikhonov weight, relative to the basis' own largest singular value, that +keeps the projection from answering a deflection it cannot represent with weights far +larger than the airfoil they correct. A plain pseudoinverse has no such bound: a kinked +deflection — a strut buckling is one — lands on the basis' weakest directions and comes +back amplified a hundredfold and alternating in sign, which is not an airfoil. The ridge +trades a small, measured under-response on deflections the basis *can* hold for a bounded +answer on the ones it cannot. + +CST is linear in its weights, so a surface displacement is a matvec against this +constant matrix — never a refit. Refitting inside a loop is not an option: the Kulfan +fit is non-unique, so the same shape fitted twice returns weight vectors differing by +more than the deformation signal, which reaches the polars as frame-to-frame jitter. + +`C(1) = 0`, so the parameterisation cannot put the trailing edge off the chord line. +Deflections must therefore be measured **against the deformed chord**, with the chord +rotation and stretch taken from the leading- and trailing-edge points; what is left is +the representable residual. +""" +struct KulfanBasis + "Chord stations the deflection is sampled on, ascending in `[0, 1]`." + x::Vector{Float64} + "`n_weights × length(x)` regularised inverse of `C(x)·B(x)`, mapping deflection to weights." + projection::Matrix{Float64} + "Number of CST weights per surface, matching the airfoil being deformed." + n_weights::Int + "Tikhonov weight the projection was built with, relative to the basis' largest singular value." + ridge::Float64 +end + +function KulfanBasis(; n_stations::Int=60, n_weights::Int=8, ridge::Real=1e-3) + n_stations > n_weights || throw(ArgumentError( + "KulfanBasis needs more stations than weights; got $n_stations and $n_weights.")) + ridge >= 0 || throw(ArgumentError("KulfanBasis ridge must not be negative.")) + theta = range(0, pi, n_stations) + x = @. (1 - cos(theta)) / 2 + shape = class_function(x) .* bernstein_basis(x, n_weights - 1) + scale = maximum(svdvals(shape))^2 + projection = iszero(ridge) ? pinv(shape) : + (shape' * shape + ridge * scale * I) \ shape' + return KulfanBasis(collect(x), projection, n_weights, Float64(ridge)) +end + +""" + deform_kulfan(basis, base, upper_deflection, lower_deflection) -> KulfanParameters + deform_kulfan(basis, base, camber) -> KulfanParameters + +Add a surface deflection to a fixed set of Kulfan parameters, both deflections sampled +on `basis.x` and normalized by chord. The three-argument form applies one camber +deflection to both surfaces, leaving the thickness distribution untouched; the +four-argument form deforms the surfaces independently, which is what a double-skin +membrane with its own upper and lower control points needs. + +Each deflection is first reduced to the part the basis can carry, see +[`chord_residual`](@ref): the straight line through its own endpoints is a chord +rotation and translation, which the basis cannot express and `pinv` answers with +runaway weights. Take that line from [`chord_line`](@ref) if the frame the deflection +was measured in has not already absorbed it. + +The leading-edge weight and the trailing-edge thickness are carried over unchanged: +they are the two shape freedoms a chord-referenced deflection cannot resolve. +""" +function deform_kulfan(basis::KulfanBasis, base::KulfanParameters, + upper_deflection::AbstractVector, + lower_deflection::AbstractVector) + length(base.upper_weights) == basis.n_weights || throw(ArgumentError( + "KulfanBasis has $(basis.n_weights) weights, airfoil has " * + "$(length(base.upper_weights)).")) + (length(upper_deflection) == length(basis.x) && + length(lower_deflection) == length(basis.x)) || throw(ArgumentError( + "Deflections must be sampled on the basis' $(length(basis.x)) stations.")) + return KulfanParameters( + base.upper_weights .+ basis.projection * chord_residual(basis, upper_deflection), + base.lower_weights .+ basis.projection * chord_residual(basis, lower_deflection), + base.leading_edge_weight, base.TE_thickness) +end + +""" + chord_residual(basis, deflection) -> Vector{Float64} + chord_line(basis, deflection) -> (offset, slope) + +The part of a deflection the CST basis can represent, and the part it cannot. `C(x)` +vanishes at both chord ends, so the basis can put neither the leading nor the trailing +edge off the chord line: a deflection that does not end at zero is asking for a chord +**rotation and translation**, not a camber change. Projecting it anyway is not merely +inexact, it is unstable — `pinv` answers an unrepresentable end displacement with +weights an order of magnitude past the ones it is correcting, and the airfoil that +comes back is not one. + +[`chord_residual`](@ref) removes the straight line through the deflection's own +endpoints and returns what is left. [`chord_line`](@ref) returns that line as +`(offset, slope)`, both over chord: `offset` moves the leading edge and `slope` is the +chord rotation `atan(slope)` [rad] the caller owes its angle of attack, when the frame +it measured the deflection in has not already absorbed it. +""" +function chord_residual(basis::KulfanBasis, deflection::AbstractVector) + offset, slope = chord_line(basis, deflection) + return deflection .- (offset .+ slope .* basis.x) +end + +function chord_line(basis::KulfanBasis, deflection::AbstractVector) + length(deflection) == length(basis.x) || throw(ArgumentError( + "Deflection must be sampled on the basis' $(length(basis.x)) stations.")) + span = basis.x[end] - basis.x[1] + abs(span) < eps() && return (deflection[1], 0.0) + slope = (deflection[end] - deflection[1]) / span + return (deflection[1] - slope * basis.x[1], slope) +end + +deform_kulfan(basis::KulfanBasis, base::KulfanParameters, + camber::AbstractVector) = deform_kulfan(basis, base, camber, camber) + +""" + control_point_deflection(basis, fractions, deflections) -> Vector{Float64} + +Resample a deflection given at arbitrary chord `fractions` onto `basis.x`, by linear +interpolation, held flat outside the sampled range. This is the generic bridge from +control points — beam nodes now, membrane nodes later — to the CST basis: the caller +only has to say where along the chord each point sits and how far off the chord line it +moved, both normalized by the local chord. + +`fractions` need not be sorted, but must not repeat a station. +""" +function control_point_deflection(basis::KulfanBasis, fractions::AbstractVector, + deflections::AbstractVector) + length(fractions) == length(deflections) || throw(ArgumentError( + "control_point_deflection: $(length(fractions)) fractions for " * + "$(length(deflections)) deflections.")) + order = sortperm(collect(float.(fractions))) + knots = collect(float.(fractions))[order] + values = collect(float.(deflections))[order] + if length(knots) < 2 + return fill(isempty(values) ? 0.0 : values[1], length(basis.x)) + end + allunique(knots) || throw(ArgumentError( + "control_point_deflection: repeated chord fraction in $knots.")) + interp = linear_interpolation(knots, values; extrapolation_bc=Flat()) + return [interp(xi) for xi in basis.x] +end diff --git a/src/airfoil_aero/kulfan.jl b/src/airfoil_aero/kulfan.jl index aef6656d..dce90044 100644 --- a/src/airfoil_aero/kulfan.jl +++ b/src/airfoil_aero/kulfan.jl @@ -6,24 +6,6 @@ This module provides functions to: - Generate airfoil coordinates from Kulfan parameters """ -""" - KulfanParameters - -Kulfan CST parameters for an airfoil. - -# Fields -- `upper_weights::Vector{Float64}`: 8 weights for upper surface -- `lower_weights::Vector{Float64}`: 8 weights for lower surface -- `leading_edge_weight::Float64`: Leading edge modification weight -- `TE_thickness::Float64`: Trailing edge thickness -""" -struct KulfanParameters - upper_weights::Vector{Float64} - lower_weights::Vector{Float64} - leading_edge_weight::Float64 - TE_thickness::Float64 -end - """ KulfanFitMethod diff --git a/src/airfoil_aero/live_polar.jl b/src/airfoil_aero/live_polar.jl new file mode 100644 index 00000000..65e37515 --- /dev/null +++ b/src/airfoil_aero/live_polar.jl @@ -0,0 +1,313 @@ +""" + LivePolarSettings(; offsets=deg2rad.(-12.0:3.0:12.0), model_size="xlarge", + weights_dir=nothing, n_crit=9.0) + +How a live polar is sampled. Every solve, each panel's deformed shape is evaluated at +its reference angle of attack plus each of `offsets`, and those values become the +panel's `SAMPLED` polar directly — there is no fit in between. + +Sampling rather than fitting is what lets a panel hold a stall. A polynomial over the +same window averages the knee into a slope and past the peak returns a lift slope of +the wrong sign, which is a divergence the solve cannot recover from; sampled values +are the polar at their own angles and reproduce whatever shape lies between them. + +The offsets set both what the polar resolves and how far the solve may move before it +reads a flat end (see [`polar_drift`](@ref)): their spacing is the resolution a stall +knee is caught at, their reach is the room the solve has. The default spacing matches +the 3° grid the offline polar tables are generated on, over a range wide enough that a +tip crossing into stall stays inside it. +""" +@with_kw struct LivePolarSettings + "Angles of attack [rad] off the reference angle, ascending and straddling zero." + offsets::Vector{Float64} = deg2rad.(collect(-12.0:3.0:12.0)) + "NeuralFoil network size." + model_size::String = "xlarge" + "Directory holding the network weights; `nothing` takes the packaged ones." + weights_dir::Union{Nothing, String} = nothing + "Critical amplification factor of the transition model." + n_crit::Float64 = 9.0 +end + +""" + LivePolars(base; settings=LivePolarSettings(), n_stations=60) + +Live polar source for a wing whose panels each carry one undeformed airfoil in `base`. +Holds the fixed CST basis, the per-panel reference angles and the network input scratch. +A refresh is dominated by the forward pass itself; the deformation is a matvec against +the constant basis, about half a microsecond a panel. Drive it with +[`refresh_live_polars!`](@ref). +""" +mutable struct LivePolars + "Sampling configuration." + settings::LivePolarSettings + "The fixed CST basis every panel's deflection is projected onto." + basis::KulfanBasis + "Undeformed Kulfan parameters per panel." + base::Vector{KulfanParameters} + "Deformed Kulfan parameters per panel, rewritten every refresh." + deformed::Vector{KulfanParameters} + "Reference angle [rad] the last refresh sampled about, per panel." + alpha_ref::Vector{Float64} + "One panel's sample angles [rad], rebuilt per panel inside a refresh." + knots::Vector{Float64} + "`25 × (n_panels · n_samples)` network input scratch." + inputs::Matrix{Float32} + "`25 × n_panels` network input scratch for the surface-pressure pass." + pressure_inputs::Matrix{Float32} +end + +function LivePolars(base::AbstractVector{KulfanParameters}; + settings::LivePolarSettings=LivePolarSettings(), + n_stations::Int=60) + offsets = settings.offsets + length(offsets) >= 2 || throw(ArgumentError( + "LivePolars needs at least two sample offsets; got $(length(offsets)).")) + issorted(offsets) || throw(ArgumentError( + "LivePolars needs ascending sample offsets.")) + first(offsets) <= 0 <= last(offsets) || throw(ArgumentError( + "LivePolars sample offsets must straddle zero, so the reference angle lies " * + "inside the sampled range; got $(rad2deg.(extrema(offsets))) deg.")) + n_panels = length(base) + return LivePolars(settings, KulfanBasis(; n_stations, + n_weights=length(base[1].upper_weights)), + collect(base), collect(base), zeros(n_panels), + zeros(length(offsets)), + zeros(Float32, 25, n_panels * length(offsets)), + zeros(Float32, 25, n_panels)) +end + +""" + panel_kulfan_parameters(panels; delta=0.0) -> Vector{KulfanParameters} + +Fit the undeformed Kulfan parameters of every panel from the surface contour its +`section_aero` carries, the starting point [`LivePolars`](@ref) deforms. + +The contour is [`shrink_wrap`](@ref)ped first, the same route the offline polar +generator takes: fitting a raw slice directly returns weights that oscillate by an +order of magnitude more than the deformation ever will, and neighbouring panels then +disagree enough to cost the VSM solve its convergence. The wrap costs ~15 ms a panel, +which is why this is done once at build time and never inside the loop. +""" +function panel_kulfan_parameters(panels; delta=0.0) + return map(panels) do panel + panel.section_aero === nothing && throw(ArgumentError( + "Live polars need a section_aero contour on every panel.")) + x, y, _, _ = section_surface(panel.section_aero, 0.0, delta) + xw, yw = shrink_wrap(collect(float.(x)), collect(float.(y)), ShrinkWrap()) + return fit_kulfan_parameters(xw, yw, LeastSquaresFit()) + end +end + +""" + polar_drift(live, alpha) -> Float64 + +How far the largest panel angle of attack has drifted off the reference angle its polar +was sampled about, as a fraction of the reach the samples give it. Above `1` a panel is +being evaluated past the last sample, where the polar is held flat and says nothing +about how the panel is really behaving, so this is the diagnostic a caller's solve loop +reports or refreshes on. Asymmetric offsets are measured by their shorter side. +""" +function polar_drift(live::LivePolars, alpha::AbstractVector) + length(alpha) == length(live.alpha_ref) || throw(ArgumentError( + "polar_drift: $(length(alpha)) angles for $(length(live.alpha_ref)) panels.")) + offsets = live.settings.offsets + reach = min(-first(offsets), last(offsets)) + return maximum(abs.(alpha .- live.alpha_ref)) / reach +end + +""" + deform_live_shapes!(live::LivePolars, deflection) -> Vector{KulfanParameters} + +Deform every base airfoil by its camber increment and store the result in +`live.deformed`, which is returned. `nothing` leaves the base shapes in place. +The deformation is an analytic perturbation of one fixed weight vector, so it +never refits and never inherits the non-uniqueness of a fit. +""" +function deform_live_shapes!(live::LivePolars, deflection) + for i in eachindex(live.base) + live.deformed[i] = isnothing(deflection) ? live.base[i] : + deform_kulfan(live.basis, live.base[i], deflection[i]) + end + return live.deformed +end + +""" + apply_live_shapes!(live::LivePolars, panels; deflection=nothing) + +Write each panel's deformed airfoil into `live_shape` without touching its polar, +so the shape a panel reports is the one its current deformation gives. Replaying a +log re-derives the drawn airfoil this way: the frame is never solved, so its polars +would say nothing, and the network pass they cost is the whole price of a refresh. +""" +function apply_live_shapes!(live::LivePolars, panels; deflection=nothing) + length(panels) == length(live.base) || throw(ArgumentError( + "apply_live_shapes!: $(length(panels)) panels for $(length(live.base)) " * + "base airfoils.")) + deform_live_shapes!(live, deflection) + for (panel, shape) in zip(panels, live.deformed) + panel.live_shape = shape + end + return nothing +end + +""" + refresh_live_polars!(live, panels, alpha_ref, reynolds; deflection=nothing) + -> Float64 + +Regenerate every panel's polar from its current shape and write it in as a `SAMPLED` +polar (see [`set_sampled_polar!`](@ref VortexStepMethod.set_sampled_polar!)). Per panel: +deform the base airfoil by `deflection` (a chord-normalized deflection on `live.basis.x`, +or `nothing` to keep the base shape), evaluate NeuralFoil at `alpha_ref .+ offsets`, and +hand those values straight to the panel. + +The write is in place — same knot count, same vectors — so a refresh every solve costs +the forward pass and nothing else. + +Each panel keeps the deformed shape it was evaluated at as its `live_shape`, so a plot +draws the airfoil the network actually saw. + +`alpha_ref` [rad] and `reynolds` are per panel or scalar. Every panel and sample goes +through the network in one forward pass. Returns the lowest analysis confidence over +the batch, a value near zero meaning the deformed shape has left the region the network +was trained on. + +Call it after the mesh has been rebuilt from the structure and before the solve: a mesh +rebuild re-seeds each panel's aero from its section, which would drop the polar. +""" +function refresh_live_polars!(live::LivePolars, panels, alpha_ref, reynolds; + deflection=nothing) + n_panels = length(live.base) + length(panels) == n_panels || throw(ArgumentError( + "refresh_live_polars!: $(length(panels)) panels for $n_panels base airfoils.")) + per_panel(v) = v isa Number ? fill(float(v), n_panels) : collect(float.(v)) + alpha_vec, re_vec = per_panel(alpha_ref), per_panel(reynolds) + live.alpha_ref .= alpha_vec + + offsets = live.settings.offsets + n_samples = length(offsets) + deform_live_shapes!(live, deflection) + for i in 1:n_panels + for k in 1:n_samples + fill_case_input!(live.inputs, (i - 1) * n_samples + k, live.deformed[i], + rad2deg(alpha_vec[i] + offsets[k]), re_vec[i], + live.settings.n_crit, 1.0, 1.0) + end + end + + model = load_neuralfoil_model(live.settings.model_size; + weights_dir=live.settings.weights_dir) + cl, cd, cm, confidence = decode_coefficients(fused_output(live.inputs, model)) + + for i in 1:n_panels + samples = ((i - 1) * n_samples + 1):(i * n_samples) + live.knots .= alpha_vec[i] .+ offsets + set_sampled_polar!(panels[i], live.knots, view(cl, samples), + view(cd, samples), view(cm, samples); + shape=live.deformed[i]) + end + return minimum(confidence) +end + +""" + refresh_live_pressure!(cp, live, contour_x, leading_edge, alpha, reynolds) + -> Vector{Vector{Float64}} + +Fill `cp[i]` with the surface pressure of panel `i`'s current deformed shape at +`alpha[i]`, resampled onto the contour nodes `contour_x[i]` (see +[`contour_pressure`](@ref)). `cp` is written in place and returned. + +This is the pressure half of a live polar. [`refresh_live_polars!`](@ref) makes the +panel forces follow the deformed shape; without this the pattern that spreads those +forces over the structure would still come from the undeformed section, so a +deformation would change how hard a panel pulls but not where it pulls. + +One batched forward pass over all panels, at the converged angle of attack rather +than at the sampled ones — a panel's `Cp` is wanted at exactly one angle, and +evaluating there is both cheaper than storing the samples and exact. Call it after +the solve has converged, with `contour_x` and `leading_edge` from the panel contours +the traction pattern is indexed on. Deform the shapes first, which +[`refresh_live_polars!`](@ref) already did for this solve. +""" +function refresh_live_pressure!(cp, live::LivePolars, contour_x, leading_edge, + alpha, reynolds) + n_panels = length(live.base) + length(cp) == length(contour_x) == length(leading_edge) == n_panels || + throw(ArgumentError("refresh_live_pressure!: $(length(cp)) pressure and " * + "$(length(contour_x)) contour entries for $n_panels panels.")) + per_panel(v) = v isa Number ? fill(float(v), n_panels) : collect(float.(v)) + alpha_vec, re_vec = per_panel(alpha), per_panel(reynolds) + for i in 1:n_panels + fill_case_input!(live.pressure_inputs, i, live.deformed[i], + rad2deg(alpha_vec[i]), re_vec[i], live.settings.n_crit, + 1.0, 1.0) + end + model = load_neuralfoil_model(live.settings.model_size; + weights_dir=live.settings.weights_dir) + station_x, cp_upper, cp_lower = decode_surface_pressure( + fused_output(live.pressure_inputs, model)) + for i in 1:n_panels + cp[i] .= contour_pressure(station_x, view(cp_upper, :, i), + view(cp_lower, :, i), contour_x[i], + leading_edge[i]) + end + return cp +end + +""" + live_surface_friction!(cf, contour_x, reynolds) -> Vector{Vector{Float64}} + +Fill `cf[i]` with the skin friction of panel `i`'s contour nodes at `reynolds[i]`, +by the same flat-plate closure ([`flat_plate_cf`](@ref)) the offline tables carry — +NeuralFoil does not predict skin friction. It depends on chord fraction and Reynolds +only, not on the shape or the angle of attack, so the live value differs from the +tabulated one purely by being at the panel's own flight Reynolds instead of the one +the tables were generated at. +""" +function live_surface_friction!(cf, contour_x, reynolds) + length(cf) == length(contour_x) || throw(ArgumentError( + "live_surface_friction!: $(length(cf)) friction and " * + "$(length(contour_x)) contour entries.")) + re_vec = reynolds isa Number ? fill(float(reynolds), length(cf)) : + collect(float.(reynolds)) + for i in eachindex(cf) + cf[i] .= (flat_plate_cf(clamp(x, 0.0, 1.0), re_vec[i]) + for x in contour_x[i]) + end + return cf +end + +""" + contour_shape_matrix(contour_x, n_weights) -> Matrix{Float64} + +The CST shape matrix `C(x)·B(x)` at chord fractions `contour_x`, mapping a change in +Kulfan weights to the normal offset it produces there. Build it once per contour: the +chord fractions of a panel's contour nodes never move, only the weights do. +""" +function contour_shape_matrix(contour_x, n_weights) + x = clamp.(collect(float.(contour_x)), 0.0, 1.0) + return class_function(x) .* bernstein_basis(x, n_weights - 1) +end + +""" + live_shape_offset!(offset, live::LivePolars, shape) -> Vector{Vector{Float64}} + +Fill `offset[i]` with the normal offset, over chord, that panel `i`'s current +deformation adds to its contour, given the panel's `shape` matrix from +[`contour_shape_matrix`](@ref). Written in place and returned. + +A deflection deforms both surfaces by the same camber increment, so one offset serves +the upper and lower halves of a contour alike and the surface split +[`contour_pressure`](@ref) needs does not arise here. Added to a reference contour it +gives the surface the network was actually evaluated on, which is what a traction +pattern has to be draped over for its normals and segment areas to mean anything. +""" +function live_shape_offset!(offset, live::LivePolars, shape) + length(offset) == length(shape) == length(live.base) || throw(ArgumentError( + "live_shape_offset!: $(length(offset)) offset and $(length(shape)) shape " * + "entries for $(length(live.base)) panels.")) + for i in eachindex(live.base) + mul!(offset[i], shape[i], + live.deformed[i].upper_weights .- live.base[i].upper_weights) + end + return offset +end diff --git a/src/airfoil_aero/neuralfoil.jl b/src/airfoil_aero/neuralfoil.jl index aea9d062..afd9928e 100644 --- a/src/airfoil_aero/neuralfoil.jl +++ b/src/airfoil_aero/neuralfoil.jl @@ -165,11 +165,40 @@ function squared_mahalanobis_distance(x::AbstractMatrix, model::NeuralFoilModel) return result end +""" + fill_case_input!(x, case, params::KulfanParameters, alpha_deg, Re, n_crit, + xtr_upper, xtr_lower) + +Write one NeuralFoil input column: rows 1–18 the Kulfan shape, rows 19–25 the flow +condition (`alpha_deg` in degrees). The single point where the network's input layout +is defined, shared by every `prepare_inputs` method. +""" +function fill_case_input!(x::AbstractMatrix, case::Int, params::KulfanParameters, + alpha_deg, Re, n_crit, xtr_upper, xtr_lower) + for i in 1:8 + x[i, case] = params.upper_weights[i] + x[8 + i, case] = params.lower_weights[i] + end + x[17, case] = params.leading_edge_weight + x[18, case] = params.TE_thickness * 50 # Scale factor from NeuralFoil + x[19, case] = sind(2 * alpha_deg) + x[20, case] = cosd(alpha_deg) + x[21, case] = 1 - cosd(alpha_deg)^2 + x[22, case] = (log(Re) - 12.5) / 3.5 + x[23, case] = (n_crit - 9) / 4.5 + x[24, case] = xtr_upper + x[25, case] = xtr_lower + return nothing +end + """ prepare_inputs(params::KulfanParameters, alpha, Re; n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) + prepare_inputs(params::AbstractVector{KulfanParameters}, alpha, Re; kwargs...) -Prepare neural network inputs from Kulfan parameters and flow conditions. +Prepare neural network inputs from Kulfan parameters and flow conditions. The vector +form takes one shape per case, so a whole wing's panels go through the network in a +single forward pass. # Arguments - `params`: Kulfan CST parameters @@ -184,39 +213,27 @@ Prepare neural network inputs from Kulfan parameters and flow conditions. """ function prepare_inputs(params::KulfanParameters, alpha, Re; n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) - # Ensure vectors alpha_vec = alpha isa Number ? [alpha] : collect(alpha) - Re_vec = Re isa Number ? fill(Re, length(alpha_vec)) : collect(Re) n_cases = length(alpha_vec) + return prepare_inputs(fill(params, n_cases), alpha_vec, Re; + n_crit, xtr_upper, xtr_lower) +end - # Broadcast other parameters - n_crit_vec = n_crit isa Number ? fill(n_crit, n_cases) : collect(n_crit) - xtr_upper_vec = xtr_upper isa Number ? fill(xtr_upper, n_cases) : collect(xtr_upper) - xtr_lower_vec = xtr_lower isa Number ? fill(xtr_lower, n_cases) : collect(xtr_lower) +function prepare_inputs(params::AbstractVector{KulfanParameters}, alpha, Re; + n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) + alpha_vec = alpha isa Number ? fill(alpha, length(params)) : collect(alpha) + n_cases = length(alpha_vec) + length(params) == n_cases || throw(ArgumentError( + "prepare_inputs: $(length(params)) shapes for $n_cases angles of attack.")) + per_case(v) = v isa Number ? fill(v, n_cases) : collect(v) + Re_vec, n_crit_vec = per_case(Re), per_case(n_crit) + xtr_upper_vec, xtr_lower_vec = per_case(xtr_upper), per_case(xtr_lower) - # Build input matrix (25 x n_cases) x = zeros(Float32, 25, n_cases) - - # Kulfan parameters (same for all cases) - for i in 1:8 - x[i, :] .= params.upper_weights[i] - x[8+i, :] .= params.lower_weights[i] - end - x[17, :] .= params.leading_edge_weight - x[18, :] .= params.TE_thickness * 50 # Scale factor from NeuralFoil - - # Flow conditions (vary per case) for i in 1:n_cases - a = alpha_vec[i] - x[19, i] = sind(2 * a) - x[20, i] = cosd(a) - x[21, i] = 1 - cosd(a)^2 - x[22, i] = (log(Re_vec[i]) - 12.5) / 3.5 - x[23, i] = (n_crit_vec[i] - 9) / 4.5 - x[24, i] = xtr_upper_vec[i] - x[25, i] = xtr_lower_vec[i] + fill_case_input!(x, i, params[i], alpha_vec[i], Re_vec[i], n_crit_vec[i], + xtr_upper_vec[i], xtr_lower_vec[i]) end - return x end @@ -323,15 +340,22 @@ function neuralfoil_aero(params::KulfanParameters, alpha, Re; n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) y = neuralfoil_fused_output(params, alpha, Re; model_size, weights_dir, n_crit, xtr_upper, xtr_lower) - analysis_confidence = sigmoid.(y[1, :]) - CL = y[2, :] ./ 2 - CD = clamp.(exp.((y[3, :] .- 2) .* 2), 0.0, 1.0) - CM = y[4, :] ./ 20 - + cl, cd, cm, confidence = decode_coefficients(y) alpha_vec = alpha isa Number ? [Float64(alpha)] : Float64.(collect(alpha)) + return NeuralFoilResult(alpha_vec, cl, cd, cm, confidence) +end + +""" + decode_coefficients(y) -> (cl, cd, cm, confidence) - return NeuralFoilResult(alpha_vec, Vector{Float64}(CL), Vector{Float64}(CD), - Vector{Float64}(CM), Vector{Float64}(analysis_confidence)) +Turn a fused network output matrix into the integrated coefficients per case, undoing +NeuralFoil's output scaling. The single place that scaling is written down. +""" +function decode_coefficients(y::AbstractMatrix) + return (Vector{Float64}(y[2, :] ./ 2), + Vector{Float64}(clamp.(exp.((y[3, :] .- 2) .* 2), 0.0, 1.0)), + Vector{Float64}(y[4, :] ./ 20), + Vector{Float64}(sigmoid.(y[1, :]))) end """ @@ -343,11 +367,21 @@ top/bottom-flipped case, flip its outputs back, and average. Returns the fused output matrix (`n_outputs × n_cases`), with the Mahalanobis penalty already applied to the confidence logit (row 1). """ -function neuralfoil_fused_output(params::KulfanParameters, alpha, Re; +function neuralfoil_fused_output(params, alpha, Re; model_size::String="xlarge", weights_dir=nothing, n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) - model = load_neuralfoil_model(model_size; weights_dir) x = prepare_inputs(params, alpha, Re; n_crit, xtr_upper, xtr_lower) + return fused_output(x, load_neuralfoil_model(model_size; weights_dir)) +end + +""" + fused_output(x, model) -> Matrix + +Symmetry-fused forward pass over a prepared input matrix, see +[`neuralfoil_fused_output`](@ref). Takes the inputs already built so a caller that +assembles its own batch does not go back through [`prepare_inputs`](@ref). +""" +function fused_output(x::AbstractMatrix, model::NeuralFoilModel) y = nn_forward(x, model) y[1, :] .-= squared_mahalanobis_distance(x, model) ./ (2 * model.n_inputs) @@ -373,18 +407,46 @@ function neuralfoil_section(params::KulfanParameters, alpha, Re; n_crit=9.0, xtr_upper=1.0, xtr_lower=1.0) y = neuralfoil_fused_output(params, alpha, Re; model_size, weights_dir, n_crit, xtr_upper, xtr_lower) - N = (size(y, 1) - 6) ÷ 6 - upper_ue = y[(7 + 2N):(6 + 3N), :] - lower_ue = y[(7 + 5N):(6 + 6N), :] alpha_vec = alpha isa Number ? [Float64(alpha)] : Float64.(collect(alpha)) + station_x, cp_upper, cp_lower = decode_surface_pressure(y) return (; alpha = alpha_vec, cl = Vector{Float64}(y[2, :] ./ 2), cd = Vector{Float64}(clamp.(exp.((y[3, :] .- 2) .* 2), 0.0, 1.0)), cm = Vector{Float64}(y[4, :] ./ 20), confidence = Vector{Float64}(sigmoid.(y[1, :])), - x = compute_optimal_x_points(N), - cp_upper = Matrix{Float64}(1 .- upper_ue .^ 2), - cp_lower = Matrix{Float64}(1 .- lower_ue .^ 2)) + x = station_x, cp_upper, cp_lower) +end + +""" + decode_surface_pressure(y) -> (station_x, cp_upper, cp_lower) + +Surface pressure of a fused output matrix, reconstructed from the predicted +edge-velocity ratios as `Cp = 1 - (ue/vinf)^2`. `station_x` is NeuralFoil's own +`N` fixed stations, and each `Cp` matrix is `N × n_cases`. +""" +function decode_surface_pressure(y::AbstractMatrix) + N = (size(y, 1) - 6) ÷ 6 + upper_ue = y[(7 + 2N):(6 + 3N), :] + lower_ue = y[(7 + 5N):(6 + 6N), :] + return (compute_optimal_x_points(N), + Matrix{Float64}(1 .- upper_ue .^ 2), + Matrix{Float64}(1 .- lower_ue .^ 2)) +end + +""" + contour_pressure(station_x, cp_upper, cp_lower, contour_x, leading_edge) + -> Vector{Float64} + +`Cp` at every node of a closed contour, taken from the two surface distributions +NeuralFoil predicts at `station_x`. Nodes up to `leading_edge` read the upper +surface and the rest the lower, which is the Selig ordering the contour is stored +in — trailing edge, over the top to the nose, back along the bottom. +""" +function contour_pressure(station_x, cp_upper, cp_lower, contour_x, leading_edge) + upper = linear_interpolation(station_x, cp_upper; extrapolation_bc=Line()) + lower = linear_interpolation(station_x, cp_lower; extrapolation_bc=Line()) + return [(k <= leading_edge ? upper : lower)(clamp(contour_x[k], 0.0, 1.0)) + for k in eachindex(contour_x)] end """ diff --git a/src/body_aerodynamics.jl b/src/body_aerodynamics.jl index a1f3cdb9..a97b007e 100644 --- a/src/body_aerodynamics.jl +++ b/src/body_aerodynamics.jl @@ -224,6 +224,13 @@ function calculate_stall_angle_list!(stall_angles::AbstractVector, # Default stall angle if none found panel_stall = stall_angle_if_none_detected + # A live polar only spans the window it was sampled over, so neither its + # curvature nor its flat ends may be read as a stall peak outside that. + if panel.aero_model == SAMPLED + stall_angles[idx] = panel_stall + continue + end + # Start with minimum cl cl_old = cl_initial diff --git a/src/panel.jl b/src/panel.jl index 83edc848..0266eab0 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -36,6 +36,10 @@ Represents a panel in a vortex step method simulation. All points and vectors ar ): Panel filaments, see: [BoundFilament](@ref) - `delta`::T=0: flap trailing-edge deflection [rad] - `crease_frac`::T=0: chordwise flap-hinge fraction (0–1); 0 disables the plate kink +- `alpha_ref`::Float64=0: reference angle [rad] of the `SAMPLED` polar +- `alpha_window`::Float64=0: half width [rad] that polar reaches; 0 = unbounded +- `alpha_knots`::Vector{Float64}=Float64[]: ascending angles [rad] a `SAMPLED` polar holds values at +- `live_shape`::Union{Nothing, KulfanParameters}=nothing: the deformed airfoil the polar was generated from """ @with_kw mutable struct Panel{T, CL, CD, CM, SA} TE_point_1::MVector{3, T} = zeros(MVector{3, T}) @@ -70,6 +74,10 @@ Represents a panel in a vortex step method simulation. All points and vectors ar ) delta::T = zero(T) crease_frac::T = zero(T) + alpha_ref::Float64 = 0.0 + alpha_window::Float64 = 0.0 + alpha_knots::Vector{Float64} = Float64[] + live_shape::Union{Nothing, KulfanParameters} = nothing end """ @@ -224,6 +232,75 @@ function init_aero!(panel::Panel, section_1::Section, section_2::Section; return nothing end +""" + sampled_value(knots, values, alpha) + +A `SAMPLED` polar's coefficient at `alpha` [rad]: linear between the knots it was +sampled on, and the end value beyond either end. Flat ends are the point — the +polar of a stalled section has a knee that no polynomial holds, and a coefficient +carried on past the samples must stay bounded rather than follow a slope out. + +`knots` is ascending; a single knot is a constant polar. +""" +@inline function sampled_value(knots, values, alpha) + n = length(knots) + n == 0 && throw(ArgumentError("A SAMPLED polar has no knots.")) + alpha <= knots[1] && return values[1] + alpha >= knots[n] && return values[n] + hi = 2 + @inbounds while alpha > knots[hi] + hi += 1 + end + @inbounds begin + span = knots[hi] - knots[hi - 1] + blend = span > 0 ? (alpha - knots[hi - 1]) / span : zero(alpha) + return values[hi - 1] + blend * (values[hi] - values[hi - 1]) + end +end + +""" + set_sampled_polar!(panel, alphas, cl, cd, cm; shape=nothing) + +Overwrite a panel's local polar with values sampled at `alphas` [rad], ascending. +The panel's aero model is set to `SAMPLED`. + +Written in place: the knots and the three value vectors reuse the panel's own +`alpha_knots` and `cl_coeffs`/`cd_coeffs`/`cm_coeffs` storage whenever the sample +count is unchanged, which is what lets a live polar source refresh every solve +without allocating. The vectors hold sampled values here rather than polynomial +coefficients — the aero model is what says which. + +The samples are the polar, so a stall knee between two of them is represented rather +than smoothed into a slope — which is the whole reason a live source samples instead +of fitting. + +`shape` is the [`KulfanParameters`](@ref) the values were generated from, stored on +the panel as `live_shape` so a panel's polar and the shape behind it are set +together and cannot drift apart. +""" +function set_sampled_polar!(panel::Panel, alphas, cl, cd, cm; shape=nothing) + length(alphas) == length(cl) == length(cd) == length(cm) || + throw(ArgumentError("A SAMPLED polar needs one value per angle; got " * + "$(length(alphas)) angles and $(length(cl))/$(length(cd))/" * + "$(length(cm)) values.")) + issorted(alphas) || + throw(ArgumentError("A SAMPLED polar needs ascending angles.")) + panel.aero_model = SAMPLED + panel.alpha_ref = alphas[(length(alphas) + 1) ÷ 2] + panel.alpha_window = maximum(abs, alphas .- panel.alpha_ref) + panel.live_shape = shape + for (dst_sym, src) in ((:alpha_knots, alphas), (:cl_coeffs, cl), + (:cd_coeffs, cd), (:cm_coeffs, cm)) + dst = getfield(panel, dst_sym) + if length(dst) == length(src) + dst .= src + else + setfield!(panel, dst_sym, collect(Float64, src)) + end + end + return nothing +end + """ reinit!(panel, section_1, section_2, aero_center, control_point, bound_point_1, bound_point_2, x_airf, y_airf, z_airf, delta, vec; kwargs...) @@ -333,6 +410,8 @@ function calculate_cl(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} cl = 2 * cos(alpha) * sin(alpha)^2 end return R(cl) + elseif panel.aero_model == SAMPLED + return R(sampled_value(panel.alpha_knots, panel.cl_coeffs, alpha)) elseif panel.aero_model == INVISCID return R(2π * alpha) end @@ -361,6 +440,9 @@ function calculate_cd(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} return R(2 * sin(alpha)^3) end return R(evalpoly(rad2deg(alpha), panel.cd_coeffs)) + elseif panel.aero_model == SAMPLED + return R(max(zero(R), + sampled_value(panel.alpha_knots, panel.cd_coeffs, alpha))) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cd_interp = panel.cd_interp cd_interp === nothing && @@ -387,6 +469,8 @@ function calculate_cm(panel::Panel{Tp}, alpha::Ta, delta::Td) where {Tp, Ta, Td} isnan(alpha) && return R(NaN) if panel.aero_model == POLY return R(evalpoly(rad2deg(alpha), panel.cm_coeffs)) + elseif panel.aero_model == SAMPLED + return R(sampled_value(panel.alpha_knots, panel.cm_coeffs, alpha)) elseif panel.aero_model in (POLAR_VECTORS, POLAR_MATRICES) cm_interp = panel.cm_interp cm_interp === nothing && diff --git a/src/section_aero.jl b/src/section_aero.jl index 3e219101..30dee3fd 100644 --- a/src/section_aero.jl +++ b/src/section_aero.jl @@ -1,3 +1,27 @@ +""" + KulfanParameters + +Kulfan CST parameters for an airfoil: the weights of the class-shape transformation +each surface is built from, the shared leading-edge modification weight and the +trailing-edge thickness. Lives here rather than in `AirfoilAero` because a +[`Panel`](@ref) carries the shape it is currently flying (`live_shape`), which a +plot or a traction pattern reads without knowing how the shape was produced. Fit one +with `AirfoilAero.fit_kulfan_parameters`, deform one with `AirfoilAero.deform_kulfan` +and turn one into coordinates with `AirfoilAero.kulfan_to_coordinates`. + +# Fields +- `upper_weights::Vector{Float64}`: weights for upper surface +- `lower_weights::Vector{Float64}`: weights for lower surface +- `leading_edge_weight::Float64`: Leading edge modification weight +- `TE_thickness::Float64`: Trailing edge thickness +""" +struct KulfanParameters + upper_weights::Vector{Float64} + lower_weights::Vector{Float64} + leading_edge_weight::Float64 + TE_thickness::Float64 +end + """ SectionAero diff --git a/test/airfoil_aero/test_live_polar.jl b/test/airfoil_aero/test_live_polar.jl new file mode 100644 index 00000000..223b5c7d --- /dev/null +++ b/test/airfoil_aero/test_live_polar.jl @@ -0,0 +1,261 @@ +using Test +using LinearAlgebra +using Statistics +using VortexStepMethod +using VortexStepMethod.AirfoilAero +using VortexStepMethod: Panel, calculate_cl, calculate_cd, calculate_cm, + set_sampled_polar! + +@testset "Kulfan deformation" begin + basis = KulfanBasis() + base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + camber = @. 0.05 * basis.x^2 * (1 - basis.x) + deformed = deform_kulfan(basis, base, camber) + + shape = AirfoilAero.class_function(basis.x) .* + AirfoilAero.bernstein_basis(basis.x, 7) + for (weights, base_weights) in ((deformed.upper_weights, base.upper_weights), + (deformed.lower_weights, base.lower_weights)) + residual = shape * (weights .- base_weights) .- camber + @test maximum(abs, residual) < 0.01 * maximum(abs, camber) + end + @test deformed.leading_edge_weight == base.leading_edge_weight + @test deformed.TE_thickness == base.TE_thickness + + # Deforming the surfaces apart changes thickness, not just camber. + apart = deform_kulfan(basis, base, camber, -camber) + @test apart.upper_weights ≉ apart.lower_weights .+ + (base.upper_weights .- base.lower_weights) + + @test_throws ArgumentError deform_kulfan(basis, base, camber[1:end-1]) +end + +@testset "an unrepresentable chord line is removed, not projected" begin + basis = KulfanBasis() + base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + camber = @. 0.02 * basis.x * (1 - basis.x) + + # A deflection whose ends are off zero is a chord rotation and translation. The + # basis cannot hold either, and projecting one anyway used to answer a 2.5% chord + # deflection with weights four times the ones it was correcting. + tilted = camber .+ 0.03 .* basis.x .+ 0.01 + offset, slope = chord_line(basis, tilted) + @test offset ≈ 0.01 atol = 1e-12 + @test slope ≈ 0.03 atol = 1e-12 + residual = chord_residual(basis, tilted) + @test residual ≈ camber atol = 1e-12 + @test abs(residual[1]) < 1e-12 && abs(residual[end]) < 1e-12 + + # So the tilted deflection has to give the same airfoil as the pure camber, and + # weights of the order of the camber rather than of the tilt. + pure = deform_kulfan(basis, base, camber) + tilt = deform_kulfan(basis, base, tilted) + @test tilt.upper_weights ≈ pure.upper_weights + @test maximum(abs, tilt.upper_weights .- base.upper_weights) < + 10 * maximum(abs, camber) + + # A pure chord rotation is no shape change at all. + rotation = deform_kulfan(basis, base, 0.05 .* basis.x) + @test rotation.upper_weights ≈ base.upper_weights atol = 1e-12 + @test rotation.lower_weights ≈ base.lower_weights atol = 1e-12 +end + +@testset "control point deflection" begin + basis = KulfanBasis() + sampled = control_point_deflection(basis, [0.0, 0.7, 0.8, 1.0], + [0.0, 0.03, 0.025, 0.0]) + @test sampled[1] ≈ 0.0 atol = 1e-12 + @test sampled[end] ≈ 0.0 atol = 1e-12 + @test maximum(sampled) ≈ 0.03 atol = 1e-3 + + # Unsorted input is the same curve as sorted input. + shuffled = control_point_deflection(basis, [0.8, 0.0, 1.0, 0.7], + [0.025, 0.0, 0.0, 0.03]) + @test shuffled ≈ sampled + @test_throws ArgumentError control_point_deflection(basis, [0.0, 0.0], [1.0, 2.0]) +end + +@testset "SAMPLED panel polar" begin + panel = Panel{Float64}() + alphas = deg2rad.([-6.0, -3.0, 0.0, 3.0, 6.0]) + set_sampled_polar!(panel, alphas, [0.0, 0.3, 0.6, 0.9, 0.8], + [0.03, 0.02, 0.02, 0.03, 0.06], [-0.1, -0.1, -0.1, -0.1, 0.0]) + @test panel.aero_model == SAMPLED + @test panel.alpha_ref ≈ 0.0 + @test panel.alpha_window ≈ deg2rad(6.0) + + # A sample is the polar at its own angle, and between two it interpolates. + @test calculate_cl(panel, deg2rad(3.0)) ≈ 0.9 + @test calculate_cl(panel, deg2rad(1.5)) ≈ 0.75 + @test calculate_cd(panel, deg2rad(-4.5)) ≈ 0.025 + @test calculate_cm(panel, deg2rad(4.5)) ≈ -0.05 + + # Which is what a fit cannot do: the peak is held rather than averaged away. + @test calculate_cl(panel, deg2rad(6.0)) < calculate_cl(panel, deg2rad(3.0)) + + # Past either end the last sample is held, so a solve that steps out reads a + # value the network gave rather than an extrapolation that runs away. + @test calculate_cl(panel, deg2rad(40.0)) ≈ 0.8 + @test calculate_cl(panel, deg2rad(-40.0)) ≈ 0.0 + @test calculate_cd(panel, deg2rad(-40.0)) ≈ 0.03 + + # A sampled polar has no flap axis: delta must not change the answer. + @test calculate_cl(panel, deg2rad(1.5), deg2rad(9.0)) == + calculate_cl(panel, deg2rad(1.5)) + + knots, coeffs = panel.alpha_knots, panel.cl_coeffs + set_sampled_polar!(panel, alphas .+ deg2rad(2.0), [0.2, 0.5, 0.8, 1.1, 1.0], + [0.03, 0.02, 0.02, 0.03, 0.06], zeros(5)) + @test panel.alpha_knots === knots && panel.cl_coeffs === coeffs + + @test_throws ArgumentError set_sampled_polar!(panel, alphas, [0.0], [0.0], [0.0]) + @test_throws ArgumentError set_sampled_polar!(panel, reverse(alphas), zeros(5), + zeros(5), zeros(5)) +end + +@testset "live polars" begin + base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + n_panels = 3 + panels = [Panel{Float64}() for _ in 1:n_panels] + live = LivePolars(fill(base, n_panels)) + confidence = refresh_live_polars!(live, panels, deg2rad(6.0), 3e6) + @test 0.0 < confidence <= 1.0 + @test all(p -> p.aero_model == SAMPLED, panels) + + # Every sample is the network's own answer, not a fit through it. + for (k, offset) in enumerate(live.settings.offsets) + @test calculate_cl(panels[1], deg2rad(6.0) + offset) ≈ + neuralfoil_aero(base, 6.0 + rad2deg(offset), 3e6).CL[1] rtol = 1e-6 + end + # And between them it tracks a direct sweep. Two bounds, because they are two + # different claims: away from the knee the polar is nearly straight and the + # samples pin it down, while at the knee the error is set by how sharp the + # corner is rather than by the spacing — halving the spacing there buys almost + # nothing, and the value stays bounded either way, which is what the solve needs. + sweep(range_deg) = maximum(abs(calculate_cl(panels[1], deg2rad(6.0 + d)) - + neuralfoil_aero(base, 6.0 + d, 3e6).CL[1]) + for d in range_deg) + @test sweep(range(-12, 6, 37)) < 0.02 + @test sweep(range(-12, 12, 49)) < 0.08 + + @test polar_drift(live, fill(deg2rad(6.0), n_panels)) ≈ 0.0 atol = 1e-12 + @test polar_drift(live, fill(deg2rad(18.0), n_panels)) ≈ 1.0 + + # Deforming the camber up raises lift at the same angle of attack. + flat = calculate_cl(panels[1], deg2rad(6.0)) + camber = @. 0.02 * live.basis.x * (1 - live.basis.x) + knots = panels[1].alpha_knots + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; + deflection=fill(camber, n_panels)) + @test calculate_cl(panels[1], deg2rad(6.0)) > flat + @test panels[1].alpha_knots === knots # a refresh writes in place + + # The panel keeps the object that was sampled, not a copy of it: a plot of + # live_shape has to be a picture of what the solve flew. + @test panels[1].live_shape === live.deformed[1] + @test panels[1].live_shape.upper_weights != + live.base[1].upper_weights # the deformation reached it + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6) + @test panels[1].live_shape === live.deformed[1] + @test panels[1].live_shape.upper_weights ≈ live.base[1].upper_weights + + @test_throws ArgumentError refresh_live_polars!(live, panels[1:2], 0.0, 3e6) + @test_throws ArgumentError LivePolars(fill(base, 2); + settings=LivePolarSettings(; offsets=deg2rad.([1.0, 2.0, 3.0]))) + @test_throws ArgumentError LivePolars(fill(base, 2); + settings=LivePolarSettings(; offsets=deg2rad.([4.0, -4.0]))) +end + +@testset "live surface pressure" begin + base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + n_panels = 2 + panels = [Panel{Float64}() for _ in 1:n_panels] + settings = LivePolarSettings(; model_size="large") + live = LivePolars(fill(base, n_panels); settings) + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6) + + # The contour the traction pattern is spread over: Selig order, nose at argmin. + contour_x, contour_y = kulfan_to_coordinates(base) + contour_x = collect(float.(contour_x)) + leading_edge = argmin(contour_x) + cp = [zeros(length(contour_x)) for _ in 1:n_panels] + same = cp[1] + refresh_live_pressure!(cp, live, fill(contour_x, n_panels), + fill(leading_edge, n_panels), deg2rad(6.0), 3e6) + @test cp[1] === same # written in place + @test all(isfinite, cp[1]) + @test cp[1] ≈ cp[2] # same shape, same angle + + # This is the same pressure the offline tables were generated with, so the live + # path and the table generator must not be two different answers for one shape. + solver = NeuralFoilSolver(; model_size="large", n_crit=settings.n_crit) + section = analyze_section(solver, DeformedSection(base, contour_x, + collect(float.(contour_y))), deg2rad(6.0), 3e6) + @test cp[1] ≈ section.cp rtol = 1e-8 + + # Suction over the upper surface, pressure under the lower one. + @test minimum(cp[1][1:leading_edge]) < -0.5 + loading(p) = sum(view(p, (leading_edge + 1):length(p))) / + (length(p) - leading_edge) - sum(view(p, 1:leading_edge)) / + leading_edge + @test loading(cp[1]) > 0.0 + + # Bending the camber up loads the section harder, which is the whole claim: the + # pattern moves with the deformation. Read as net loading, not as peak suction — + # a cambered nose meets the flow better, so the leading-edge peak actually + # weakens (-2.64 to -2.57) while the section as a whole carries more. + camber = @. 0.02 * live.basis.x * (1 - live.basis.x) + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; + deflection=fill(camber, n_panels)) + cambered = deepcopy(cp) + refresh_live_pressure!(cambered, live, fill(contour_x, n_panels), + fill(leading_edge, n_panels), deg2rad(6.0), 3e6) + @test loading(cambered[1]) > loading(cp[1]) + + @test_throws ArgumentError refresh_live_pressure!(cp[1:1], live, + fill(contour_x, n_panels), fill(leading_edge, n_panels), 0.0, 3e6) +end + +@testset "the contour follows the deformation" begin + base = KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + panels = [Panel{Float64}() for _ in 1:2] + live = LivePolars(fill(base, 2)) + shape = [contour_shape_matrix(live.basis.x, 8) for _ in 1:2] + offset = [zeros(length(live.basis.x)) for _ in 1:2] + same = offset[1] + + # An undeformed panel's contour is the reference one, exactly. + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6) + live_shape_offset!(offset, live, shape) + @test offset[1] === same + @test all(iszero, offset[1]) + + # A camber increment comes back as itself: this is the same offset the airfoil + # the network saw was built from, so the contour and the pressures agree. + camber = @. 0.02 * live.basis.x * (1 - live.basis.x) + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; + deflection=fill(camber, 2)) + live_shape_offset!(offset, live, shape) + @test maximum(abs, offset[1] .- camber) < 0.02 * maximum(abs, camber) + @test offset[1] ≈ offset[2] + + # A chord rotation is not a shape change, so it moves no contour node. + refresh_live_polars!(live, panels, deg2rad(6.0), 3e6; + deflection=fill(0.05 .* live.basis.x, 2)) + live_shape_offset!(offset, live, shape) + @test maximum(abs, offset[1]) < 1e-12 + + @test_throws ArgumentError live_shape_offset!(offset[1:1], live, shape) +end + +@testset "live skin friction" begin + contour_x = [1.0, 0.5, 0.0, 0.5, 1.0] + cf = [zeros(5), zeros(5)] + same = cf[1] + live_surface_friction!(cf, [contour_x, contour_x], [3e6, 1e6]) + @test cf[1] === same + @test cf[1] ≈ [AirfoilAero.flat_plate_cf(x, 3e6) for x in contour_x] + # Lower Reynolds is more friction, and it is highest at the nose. + @test all(cf[2] .> cf[1]) + @test argmax(cf[1]) == 3 +end diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl index 14807a6f..afe3c9ba 100644 --- a/test/plotting/test_plotting.jl +++ b/test/plotting/test_plotting.jl @@ -521,6 +521,29 @@ end # In-place update refreshes the registered skin observables without error. @test_nowarn Makie.plot!(body_aero) + # Live polars: the skin has to draw the panel's stored deformed shape, not the + # tabulated contour, or a deformation bug is invisible in the picture. + basis = VortexStepMethod.AirfoilAero.KulfanBasis() + base = VortexStepMethod.KulfanParameters(fill(0.15, 8), fill(-0.05, 8), 0.0, 0.0) + cambered = VortexStepMethod.AirfoilAero.deform_kulfan(basis, base, + @. 0.08 * basis.x * (1 - basis.x)) + for panel in body_aero.panels + VortexStepMethod.set_sampled_polar!(panel, deg2rad.([-4.0, 0.0, 4.0]), + [0.3, 0.5, 0.7], [0.02, 0.02, 0.03], [-0.05, -0.05, -0.05]; + shape=cambered) + end + @test body_aero.panels[1].live_shape === cambered + v_live, f_live, ribs_live = airfoil_skin_geometry(body_aero) + live_nodes = length(VortexStepMethod.AirfoilAero.kulfan_to_coordinates(cambered)[1]) + @test live_nodes != n_node # the two routes are distinguishable + @test all(rib -> length(rib) == live_nodes, ribs_live) + @test length(f_live) == 2 * (n_sections - 1) * (live_nodes - 1) + @test_nowarn Makie.plot!(Axis3(Figure()[1, 1]), body_aero; airfoils=true) + for panel in body_aero.panels + panel.live_shape = nothing + end + @test all(rib -> length(rib) == n_node, airfoil_skin_geometry(body_aero)[3]) + # A body without any surface aero yields no skin geometry, but still plots. plain_body = create_body_aero() v_plain, f_plain, ribs_plain = airfoil_skin_geometry(plain_body) diff --git a/test/runtests.jl b/test/runtests.jl index d613c2b5..08f8763c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -62,6 +62,7 @@ function include_selected_tests() should_run_test("wing_geometry/test_billowing.jl") && include("wing_geometry/test_billowing.jl") should_run_test("yaml_geometry/test_yaml_geometry.jl") && include("yaml_geometry/test_yaml_geometry.jl") should_run_test("airfoil_aero/test_airfoil_aero.jl") && include("airfoil_aero/test_airfoil_aero.jl") + should_run_test("airfoil_aero/test_live_polar.jl") && include("airfoil_aero/test_live_polar.jl") should_run_test("obj_adapter/test_obj_adapter.jl") && include("obj_adapter/test_obj_adapter.jl") should_run_test("surfplan/test_surfplan.jl") && include("surfplan/test_surfplan.jl") should_run_test("Aqua.jl") && include("Aqua.jl")