diff --git a/CHANGELOG.md b/CHANGELOG.md index 9b303475..a3d41c12 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ ### Added +- `masure_regression` airfoils `{t, eta, kappa, delta, lambda, phi}` resolve to polars: `resolve_aero_geometry(...; ml_models_dir)` evaluates the Masure Extra-Trees regression (`masure_aero`, `load_masure_model`) in pure Julia, from models converted once with `scripts/export_masure_models.py`. - `plot_geometry`, `plot_distribution`, `plot_combined_analysis`, `plot_section_polars`, `plot_airfoil_fit` and `plot_airfoils` take `show_title=true`; with `false` the title is not drawn, and still names the saved file and window. diff --git a/data/TUDELFT_V3_KITE/aero_geometry.yaml b/data/TUDELFT_V3_KITE/aero_geometry.yaml index 6918cd72..a136a2f7 100644 --- a/data/TUDELFT_V3_KITE/aero_geometry.yaml +++ b/data/TUDELFT_V3_KITE/aero_geometry.yaml @@ -75,7 +75,7 @@ wing_airfoils: # - polars: # polar_file_path: Path to polar CSV file (columns: alpha [rad], cl, cd, cm) # - masure_regression: - # t, eta, kappa, delta, lamba, phi: Regression parameters + # t, eta, kappa, delta, lambda, phi: Regression parameters # - inviscid: # no further data is required # --------------------------------------------------------------- diff --git a/data/TUDELFT_V3_KITE/aero_geometry_coarse_discretisation.yaml b/data/TUDELFT_V3_KITE/aero_geometry_coarse_discretisation.yaml index f122cfbb..a1274956 100644 --- a/data/TUDELFT_V3_KITE/aero_geometry_coarse_discretisation.yaml +++ b/data/TUDELFT_V3_KITE/aero_geometry_coarse_discretisation.yaml @@ -120,7 +120,7 @@ wing_airfoils: # --------------------------------------------------------------- # headers: # - airfoil_id: integer, unique identifier for the airfoil - # - type: one of [neuralfoil, breukels_regression, measure_regression, polars] + # - type: one of [neuralfoil, breukels_regression, masure_regression, polars] # - info_dict: dictionary with parameters depending on 'type' # # info_dict fields by type: @@ -141,7 +141,7 @@ wing_airfoils: # Dirty wind tunnel: 4–8 # - polars: # polar_file_path: Path to polar CSV file (columns: alpha [rad], cl, cd, cm) - # - measure_regression: + # - masure_regression: # t, eta, kappa, delta, lambda, phi: Regression parameters # - inviscid: # no further data is required diff --git a/docs/src/airfoil_pipeline.md b/docs/src/airfoil_pipeline.md index edc5d932..9cc66fa8 100644 --- a/docs/src/airfoil_pipeline.md +++ b/docs/src/airfoil_pipeline.md @@ -124,6 +124,10 @@ keeping its own edge positions, and a warning lists the reuse. All floats are ro precision by the single [`write_yaml`](@ref VortexStepMethod.ObjAdapter.write_yaml) writer, so generated geometry files stay diff-friendly and consistent. +## Parametric LEI airfoils: the masure regression + +A geometry YAML can also describe a leading-edge-inflatable airfoil by six shape parameters instead of a slice, as the `wing_airfoils` entry `[id, masure_regression, {t, eta, kappa, delta, lambda, phi}]`: tube diameter, chordwise and vertical position of maximum camber, trailing-edge reflex angle [deg], camber tension and leading-edge tension. [`resolve_aero_geometry`](@ref VortexStepMethod.ObjAdapter.resolve_aero_geometry) turns it into a `POLAR_VECTORS` table over the block's `alpha_range` with [`masure_aero`](@ref), which evaluates the Extra-Trees regression of K.R.G. Masure, trained on 2D RANS simulations, in pure Julia. The models exist for `reynolds` 1e6, 5e6 and 2e7. Download them from [Zenodo](https://doi.org/10.5281/zenodo.16925758), convert them once with `python scripts/export_masure_models.py ET_re*.pkl --out DIR` (needs numpy and scikit-learn), and pass `ml_models_dir=DIR`. The converted model gives the same coefficients as scikit-learn's `predict`. Keep `alpha_range` within the angles the models were trained on: outside them the trees return a constant. + ## Why this is useful - **CAD in, solver-ready model out.** No hand-drawing airfoils or manually pairing them diff --git a/docs/src/functions.md b/docs/src/functions.md index 8f934a66..a1913716 100644 --- a/docs/src/functions.md +++ b/docs/src/functions.md @@ -41,6 +41,8 @@ deform_section analyze_section analyze_sweep neuralfoil_aero +masure_aero +load_masure_model deform_kulfan chord_residual chord_line diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index d8b0c836..51018425 100644 --- a/docs/src/private_functions.md +++ b/docs/src/private_functions.md @@ -228,6 +228,15 @@ swish sigmoid ``` +### Masure regression +```@docs +MASURE_PARAMETERS +MASURE_REYNOLDS +forest_predict +one_based +scale_input +``` + ### Polars and airfoil IO ```@docs create_2d_polars diff --git a/docs/src/private_types.md b/docs/src/private_types.md index 0941c427..4072b640 100644 --- a/docs/src/private_types.md +++ b/docs/src/private_types.md @@ -30,6 +30,8 @@ KulfanBasis LivePolarSettings LivePolars NeuralFoilModel +MasureModel +ExtraTreesForest NeuralFoilResult NeuralFoilWorkspace ContourPressureScratch diff --git a/scripts/export_masure_models.py b/scripts/export_masure_models.py new file mode 100644 index 00000000..25022e56 --- /dev/null +++ b/scripts/export_masure_models.py @@ -0,0 +1,125 @@ +"""Export the masure-regression scikit-learn models to the .npz files AirfoilAero reads. + + python export_masure_models.py ET_re1e6.pkl ET_re5e6.pkl ET_re2e7.pkl --out DIR + python export_masure_models.py --fixture DIR + +The first form converts the pickles from https://doi.org/10.5281/zenodo.16925758 into +`DIR/ET_re.npz`; pass `DIR` as `ml_models_dir` to `resolve_aero_geometry`. The +second trains a small model of the same layout on synthetic data and writes it with +the inputs and `predict` outputs the Julia test compares against. Needs numpy and +scikit-learn. + +An exported file holds `input_mean` and `input_scale` (the StandardScaler) and, for +output k = 0, 1, 2 (CD, CL, CM), the nodes of all its trees concatenated: +`output{k}_roots`, `_left`, `_right`, `_feature`, `_threshold`, `_value`. Node and +feature indices are 0-based, and `_left` is -1 on a leaf. +""" + +import argparse +import pickle +from pathlib import Path + +import numpy as np +from sklearn.ensemble import ExtraTreesRegressor +from sklearn.multioutput import MultiOutputRegressor +from sklearn.pipeline import Pipeline +from sklearn.preprocessing import StandardScaler + + +def forest_arrays(forest): + """Concatenate the nodes of every tree in `forest`, children indexed globally.""" + roots, left, right, feature, threshold, value = [], [], [], [], [], [] + offset = 0 + for estimator in forest.estimators_: + tree = estimator.tree_ + if tree.n_outputs != 1: + raise ValueError(f"expected single-output trees, got {tree.n_outputs}") + is_leaf = tree.children_left == -1 + roots.append(offset) + left.append(np.where(is_leaf, -1, tree.children_left + offset)) + right.append(np.where(is_leaf, -1, tree.children_right + offset)) + feature.append(tree.feature) + threshold.append(tree.threshold) + value.append(tree.value[:, 0, 0]) + offset += tree.node_count + return { + "roots": np.asarray(roots, dtype=np.int32), + "left": np.concatenate(left).astype(np.int32), + "right": np.concatenate(right).astype(np.int32), + "feature": np.concatenate(feature).astype(np.int32), + "threshold": np.concatenate(threshold).astype(np.float64), + "value": np.concatenate(value).astype(np.float64), + } + + +def export_model(model, npz_path): + """Write a StandardScaler -> MultiOutputRegressor(ExtraTreesRegressor) pipeline.""" + steps = [step for _, step in model.steps] if isinstance(model, Pipeline) else [] + if not (len(steps) == 2 and isinstance(steps[0], StandardScaler) + and isinstance(steps[1], MultiOutputRegressor)): + raise ValueError(f"unexpected model layout: {model!r}") + scaler, multi_output = steps + if len(multi_output.estimators_) != 3: + raise ValueError(f"expected 3 outputs, got {len(multi_output.estimators_)}") + arrays = {"input_mean": scaler.mean_, "input_scale": scaler.scale_} + for k, forest in enumerate(multi_output.estimators_): + if not isinstance(forest, ExtraTreesRegressor): + raise ValueError(f"output {k} is a {type(forest).__name__}") + for name, array in forest_arrays(forest).items(): + arrays[f"output{k}_{name}"] = array + np.savez(npz_path, **arrays) + + +def threshold_rows(model, row): + """Copies of `row` moved onto each output's first root split, where float32 matters.""" + scaler, multi_output = (step for _, step in model.steps) + rows = [] + for forest in multi_output.estimators_: + tree = forest.estimators_[0].tree_ + feature, threshold = tree.feature[0], tree.threshold[0] + for step in (-1e-9, 0.0, 1e-9): + moved = row.copy() + moved[feature] = threshold * (1 + step) * scaler.scale_[feature] + \ + scaler.mean_[feature] + rows.append(moved) + return np.array(rows) + + +def write_fixture(out_dir): + """Train a small pipeline of the published layout and write it with its predictions.""" + rng = np.random.default_rng(42) + low = np.array([0.05, 0.1, 0.0, -10.0, 0.1, 0.1, -10.0]) + high = np.array([0.12, 0.6, 0.15, 5.0, 0.4, 0.9, 30.0]) + X = rng.uniform(low, high, size=(400, 7)) + alpha = np.deg2rad(X[:, 6]) + y = np.column_stack([0.02 + 0.3 * alpha**2 + X[:, 0], + 2 * np.pi * alpha + 5 * X[:, 2], + -0.1 - X[:, 2] + 0.01 * X[:, 3]]) + model = Pipeline([ + ("scale", StandardScaler()), + ("model", MultiOutputRegressor(ExtraTreesRegressor( + n_estimators=5, max_depth=6, max_features="log2", random_state=42))), + ]).fit(X, y) + out_dir.mkdir(parents=True, exist_ok=True) + export_model(model, out_dir / "ET_re1e6.npz") + X_test = np.vstack([rng.uniform(low, high, size=(20, 7)), threshold_rows(model, X[0])]) + np.savez(out_dir / "reference.npz", X=X_test, Y=model.predict(X_test)) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("pickles", nargs="*", type=Path) + parser.add_argument("--out", type=Path, default=Path(".")) + parser.add_argument("--fixture", type=Path) + args = parser.parse_args() + if args.fixture is not None: + write_fixture(args.fixture) + for pickle_path in args.pickles: + with open(pickle_path, "rb") as file: + model = pickle.load(file) + args.out.mkdir(parents=True, exist_ok=True) + export_model(model, args.out / pickle_path.with_suffix(".npz").name) + + +if __name__ == "__main__": + main() diff --git a/src/airfoil_aero/AirfoilAero.jl b/src/airfoil_aero/AirfoilAero.jl index de24e654..650066f2 100644 --- a/src/airfoil_aero/AirfoilAero.jl +++ b/src/airfoil_aero/AirfoilAero.jl @@ -16,6 +16,7 @@ include("kulfan.jl") include("deform.jl") include("shrink_wrap.jl") include("neuralfoil.jl") +include("masure_regression.jl") include("poly.jl") include("airfoil_solvers/common.jl") include("airfoil_solvers/xfoil_solver.jl") @@ -33,6 +34,7 @@ export ShrinkWrap, shrink_wrap export fit_kulfan_parameters, kulfan_to_coordinates export NeuralFoilModel, NeuralFoilResult, load_neuralfoil_model export neuralfoil_aero, neuralfoil_section +export MasureModel, load_masure_model, masure_aero export KulfanBasis, deform_kulfan, control_point_deflection export chord_residual, chord_line export LivePolarSettings, LivePolars, panel_kulfan_parameters diff --git a/src/airfoil_aero/airfoil_io.jl b/src/airfoil_aero/airfoil_io.jl index b0ee6e2a..3beb97e7 100644 --- a/src/airfoil_aero/airfoil_io.jl +++ b/src/airfoil_aero/airfoil_io.jl @@ -16,15 +16,10 @@ end """ write_polar(filepath, result::NeuralFoilResult) -Write a NeuralFoil result to a `POLAR_VECTORS` table (`alpha, Cd, Cs, Cl, Cm`), CSV or -Arrow as the suffix of `filepath` says (see [`write_node_rows`](@ref -VortexStepMethod.write_node_rows)). +Write a NeuralFoil result to a `POLAR_VECTORS` table. """ -function write_polar(filepath::String, result::NeuralFoilResult) - values = [result.CD zero(result.CD) result.CL result.CM] - return write_node_rows(filepath, deg2rad.(result.alpha), nothing, values; - columns=["Cd", "Cs", "Cl", "Cm"]) -end +write_polar(filepath::String, result::NeuralFoilResult) = + write_polar(filepath, deg2rad.(result.alpha), result.CL, result.CD, result.CM) """ generate_polar_from_coordinates(x, y, output_path; Re, alpha_range=-180:1:180, @@ -69,34 +64,43 @@ function generate_polar_from_coordinates(x::Vector, y::Vector, output_path::Stri end """ - resolve_airfoil(type, info, out_dir, id; Re, alpha_range, table_format=:csv) - -> (new_type, new_info) + resolve_airfoil(type, info, out_dir, id; Re, alpha_range, table_format=:csv, + ml_models_dir=nothing) -> (new_type, new_info) Resolve one awesIO `wing_airfoils` entry to a core-loadable form. `breukels_regression` `(t, kappa)` → `poly` coeffs (via [`lei_poly_coeffs`](@ref)); `neuralfoil` -`(dat_file_path, …)` → a `polars` table `{id}.{table_format}` (`:csv` or `:arrow`, via -[`generate_polar_from_dat`](@ref)) written under `out_dir`; `polars`/`poly`/`inviscid` -pass through. `masure_regression` is not yet supported. `info` file paths should -already be absolute. +`(dat_file_path, …)` and `masure_regression` `(t, eta, kappa, delta, lambda, phi)` → a +`polars` table `{id}.{table_format}` (`:csv` or `:arrow`) written under `out_dir`, the +latter from the model in `ml_models_dir` (see [`load_masure_model`](@ref)); +`polars`/`poly`/`inviscid` pass through. `info` file paths should already be absolute. """ function resolve_airfoil(type::AbstractString, info::AbstractDict, out_dir, id; - Re, alpha_range, table_format::Symbol=:csv) + Re, alpha_range, table_format::Symbol=:csv, + ml_models_dir=nothing) + polar = joinpath(out_dir, "$(id).$(table_format)") if type == "breukels_regression" cl, cd, cm = lei_poly_coeffs(Float64(info["t"]), Float64(info["kappa"])) return "poly", Dict{String,Any}("cl_coeffs" => cl, "cd_coeffs" => cd, "cm_coeffs" => cm) elseif type == "neuralfoil" - polar = joinpath(out_dir, "$(id).$(table_format)") solver = NeuralFoilSolver( model_size=String(get(info, "model_size", "large")), n_crit=Float64(get(info, "n_crit", 9.0))) generate_polar_from_dat(String(info["dat_file_path"]), polar; Re, alpha_range=collect(alpha_range), solver) return "polars", Dict{String,Any}("polar_file_path" => polar) + elseif type == "masure_regression" + ml_models_dir === nothing && + error("masure_regression airfoil $id needs ml_models_dir.") + absent = [name for name in MASURE_PARAMETERS if !haskey(info, name)] + isempty(absent) || error("masure_regression airfoil $id lacks " * + "$(join(absent, ", ")); it takes $(join(MASURE_PARAMETERS, ", ")).") + alpha = collect(Float64, alpha_range) + cl, cd, cm = masure_aero(load_masure_model(Re, ml_models_dir), info, alpha) + write_polar(polar, deg2rad.(alpha), cl, cd, cm) + return "polars", Dict{String,Any}("polar_file_path" => polar) elseif type in ("polars", "poly", "inviscid") return String(type), info - elseif type == "masure_regression" - error("masure_regression airfoils are not yet supported by AirfoilAero.") else error("Unknown airfoil type: $type") end diff --git a/src/airfoil_aero/airfoil_solvers/common.jl b/src/airfoil_aero/airfoil_solvers/common.jl index 3a99ed13..b7fb829e 100644 --- a/src/airfoil_aero/airfoil_solvers/common.jl +++ b/src/airfoil_aero/airfoil_solvers/common.jl @@ -31,22 +31,29 @@ struct SectionSolution cf::Vector{Float64} end +""" + write_polar(filepath, alpha, cl, cd, cm) + +Write a polar to a `POLAR_VECTORS` table (`alpha, Cd, Cs, Cl, Cm`), CSV or Arrow as the +suffix of `filepath` says. `alpha` is in radians and is written in degrees. +""" +function write_polar(filepath::String, alpha, cl, cd, cm) + return write_node_rows(filepath, alpha, nothing, [cd zero(cd) cl cm]; + columns=["Cd", "Cs", "Cl", "Cm"]) +end + """ write_polar(filepath, sols::Vector{SectionSolution}) Write a solver sweep (from any [`AbstractAirfoilSolver`](@ref)) to a `POLAR_VECTORS` -table (`alpha, Cd, Cs, Cl, Cm`; alpha in degrees), CSV or Arrow as the suffix of -`filepath` says. Non-converged angles (`NaN`) are skipped, so this works for both -NeuralFoil and XFoil sweeps. +table, skipping non-converged angles (`NaN`), so this works for both NeuralFoil and +XFoil sweeps. """ function write_polar(filepath::String, sols::Vector{SectionSolution}) converged = filter(sol -> !isnan(sol.cl), sols) - alpha = [sol.alpha for sol in converged] - cd = [sol.cd for sol in converged] - cl = [sol.cl for sol in converged] - cm = [sol.cm for sol in converged] - return write_node_rows(filepath, alpha, nothing, [cd zero(cd) cl cm]; - columns=["Cd", "Cs", "Cl", "Cm"]) + return write_polar(filepath, [sol.alpha for sol in converged], + [sol.cl for sol in converged], [sol.cd for sol in converged], + [sol.cm for sol in converged]) end """ diff --git a/src/airfoil_aero/masure_regression.jl b/src/airfoil_aero/masure_regression.jl new file mode 100644 index 00000000..40d8ee29 --- /dev/null +++ b/src/airfoil_aero/masure_regression.jl @@ -0,0 +1,117 @@ +"Airfoil parameters the masure regression takes, in its input order (alpha follows)." +const MASURE_PARAMETERS = ("t", "eta", "kappa", "delta", "lambda", "phi") + +"Reynolds numbers a masure regression model exists for, with their file suffix." +const MASURE_REYNOLDS = Dict(1.0e6 => "1e6", 5.0e6 => "5e6", 2.0e7 => "2e7") + +""" + ExtraTreesForest + +The regression trees predicting one output, nodes of all trees concatenated. A leaf +has `left == 0`. +""" +struct ExtraTreesForest + "first node of each tree" + roots::Vector{Int32} + "child taken when the feature is at or below the threshold" + left::Vector{Int32} + "child taken when the feature is above the threshold" + right::Vector{Int32} + "input feature a node splits on" + feature::Vector{Int32} + "split threshold, in scaled input units" + threshold::Vector{Float64} + "prediction of a leaf" + value::Vector{Float64} +end + +""" + MasureModel + +A masure regression model for one Reynolds number: the input standardisation and one +[`ExtraTreesForest`](@ref) per output, in the order CD, CL, CM. +""" +struct MasureModel + input_mean::Vector{Float64} + input_scale::Vector{Float64} + forests::Vector{ExtraTreesForest} +end + +const MASURE_CACHE = Dict{String, MasureModel}() + +"Index array `key` of an exported model, shifted from 0-based to 1-based." +one_based(data::AbstractDict, key::String) = data[key] .+ Int32(1) + +function ExtraTreesForest(data::AbstractDict, prefix::String) + return ExtraTreesForest(one_based(data, "$(prefix)_roots"), + one_based(data, "$(prefix)_left"), + one_based(data, "$(prefix)_right"), + one_based(data, "$(prefix)_feature"), + data["$(prefix)_threshold"], data["$(prefix)_value"]) +end + +""" + load_masure_model(Re, ml_models_dir) -> MasureModel + +Load the masure regression model for Reynolds number `Re` (1e6, 5e6 or 2e7) from +`ml_models_dir/ET_re.npz`, as written by `scripts/export_masure_models.py`. +""" +function load_masure_model(Re::Real, ml_models_dir::AbstractString) + haskey(MASURE_REYNOLDS, Re) || error("No masure regression model for Re = $Re; " * + "available: $(join(sort!(collect(keys(MASURE_REYNOLDS))), ", ")).") + path = abspath(joinpath(ml_models_dir, "ET_re$(MASURE_REYNOLDS[Re]).npz")) + haskey(MASURE_CACHE, path) && return MASURE_CACHE[path] + isfile(path) || error("Masure regression model not found at $path. Download " * + "the models from https://doi.org/10.5281/zenodo.16925758 and convert them " * + "with scripts/export_masure_models.py.") + data = npzread(path) + return MASURE_CACHE[path] = MasureModel(data["input_mean"], data["input_scale"], + [ExtraTreesForest(data, "output$k") for k in 0:2]) +end + +""" + forest_predict(forest::ExtraTreesForest, x) -> Float64 + +Mean of the leaf values the scaled input `x` reaches in each tree. +""" +function forest_predict(forest::ExtraTreesForest, x::AbstractVector{Float32}) + total = 0.0 + for root in forest.roots + node = root + while forest.left[node] != 0 + node = x[forest.feature[node]] <= forest.threshold[node] ? + forest.left[node] : forest.right[node] + end + total += forest.value[node] + end + return total / length(forest.roots) +end + +""" + masure_aero(model::MasureModel, params, alpha) -> (cl, cd, cm) + +Lift, drag and moment coefficients of the LEI airfoil described by `params` (keys +[`MASURE_PARAMETERS`](@ref), `delta` in degrees) at each angle of attack in `alpha` +[deg]. +""" +function masure_aero(model::MasureModel, params::AbstractDict, alpha::AbstractVector) + x = [scale_input(model, params[name], j) for (j, name) in enumerate(MASURE_PARAMETERS)] + push!(x, 0.0f0) + coefficients = zeros(length(alpha), length(model.forests)) + for (i, angle) in enumerate(alpha) + x[end] = scale_input(model, angle, length(x)) + for (k, forest) in enumerate(model.forests) + coefficients[i, k] = forest_predict(forest, x) + end + end + return coefficients[:, 2], coefficients[:, 1], coefficients[:, 3] +end + +""" + scale_input(model, value, j) -> Float32 + +Input `j` standardised as the model's scaler does, rounded to `Float32` as sklearn trees +compare it. +""" +scale_input(model::MasureModel, value::Real, j::Int) = + Float32((Float64(value) - model.input_mean[j]) / model.input_scale[j]) diff --git a/src/obj_adapter/obj_to_yaml.jl b/src/obj_adapter/obj_to_yaml.jl index 609e84b6..84d3589e 100644 --- a/src/obj_adapter/obj_to_yaml.jl +++ b/src/obj_adapter/obj_to_yaml.jl @@ -263,17 +263,20 @@ function obj_to_yaml(obj_path::String, output_dir::String; end """ - resolve_aero_geometry(yaml_in, out_dir; table_format=:csv, verbose=true) -> yaml_out + resolve_aero_geometry(yaml_in, out_dir; table_format=:csv, ml_models_dir=nothing, + verbose=true) -> yaml_out Read an awesIO-style geometry YAML and resolve every `wing_airfoils` entry to a core-loadable form via [`resolve_airfoil`](@ref) (`breukels_regression` → `poly`, -`neuralfoil` → `polars` table, others pass through). The generated polars go to +`neuralfoil` and `masure_regression` → `polars` table, others pass through); +`ml_models_dir` holds the masure regression models. The generated polars go to `out_dir/polars` as `table_format` (`:csv` or `:arrow`), the resolved YAML to `out_dir/geometry.yaml`. `wing_sections` (incl. any `VUP` up-vectors) pass through unchanged. Load the result with `Wing(yaml_out)`. """ function resolve_aero_geometry(yaml_in::String, out_dir::String; - table_format::Symbol=:csv, verbose=true) + table_format::Symbol=:csv, ml_models_dir=nothing, + verbose=true) mkpath(out_dir) polar_dir = joinpath(out_dir, "polars") mkpath(polar_dir) @@ -296,7 +299,7 @@ function resolve_aero_geometry(yaml_in::String, out_dir::String; end new_type, new_info = resolve_airfoil(String(row[ti]), info, polar_dir, row[idi]; Re, alpha_range, - table_format) + table_format, ml_models_dir) row[ti] = new_type row[ii] = new_info end diff --git a/test/airfoil_aero/data/masure/ET_re1e6.npz b/test/airfoil_aero/data/masure/ET_re1e6.npz new file mode 100644 index 00000000..1dbcadb5 Binary files /dev/null and b/test/airfoil_aero/data/masure/ET_re1e6.npz differ diff --git a/test/airfoil_aero/data/masure/reference.npz b/test/airfoil_aero/data/masure/reference.npz new file mode 100644 index 00000000..e4b03e6e Binary files /dev/null and b/test/airfoil_aero/data/masure/reference.npz differ diff --git a/test/airfoil_aero/test_masure_regression.jl b/test/airfoil_aero/test_masure_regression.jl new file mode 100644 index 00000000..da1969e0 --- /dev/null +++ b/test/airfoil_aero/test_masure_regression.jl @@ -0,0 +1,62 @@ +using Test +import YAML +using VortexStepMethod +using VortexStepMethod.ObjAdapter: resolve_aero_geometry +using VortexStepMethod.AirfoilAero: MASURE_PARAMETERS, load_masure_model, masure_aero, + resolve_airfoil, write_yaml, npzread + +fixture_dir = joinpath(@__DIR__, "data", "masure") + +@testset "masure regression" begin + @testset "Extra-Trees evaluation matches sklearn predict, on-split rows too" begin + reference = npzread(joinpath(fixture_dir, "reference.npz")) + model = load_masure_model(1e6, fixture_dir) + for (row, expected) in zip(eachrow(reference["X"]), eachrow(reference["Y"])) + params = Dict(zip(MASURE_PARAMETERS, row[1:6])) + cl, cd, cm = masure_aero(model, params, [row[7]]) + @test [cd[1], cl[1], cm[1]] == expected + end + end + + @testset "load_masure_model rejects an untrained Re and a missing file" begin + @test_throws "No masure regression model for Re = 3.0e6" load_masure_model( + 3e6, fixture_dir) + @test_throws "Masure regression model not found" load_masure_model( + 5e6, fixture_dir) + end + + @testset "resolve_airfoil names the airfoil and the missing parameters" begin + info = Dict("t" => 0.08, "eta" => 0.2, "kappa" => 0.09, "delta" => -2.0, + "lamba" => 0.2, "phi" => 0.6) + @test_throws "masure_regression airfoil 3 lacks lambda" resolve_airfoil( + "masure_regression", info, mktempdir(), 3; Re=1e6, alpha_range=[0.0], + ml_models_dir=fixture_dir) + end + + @testset "resolve_aero_geometry turns masure_regression into a loadable polar" begin + dir = mktempdir() + yaml_in = joinpath(dir, "awesio.yaml") + params = Dict("t" => 0.08, "eta" => 0.2, "kappa" => 0.09, "delta" => -2.0, + "lambda" => 0.2, "phi" => 0.6) + write_yaml(yaml_in, Dict( + "wing_sections" => Dict( + "headers" => ["airfoil_id", "LE_x", "LE_y", "LE_z", "TE_x", "TE_y", "TE_z"], + "data" => [Any[1, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0], + Any[1, 0.0, -1.0, 0.0, 1.0, -1.0, 0.0]]), + "wing_airfoils" => Dict( + "alpha_range" => [-4, 4, 2], "reynolds" => 1e6, + "headers" => ["airfoil_id", "type", "info_dict"], + "data" => [Any[1, "masure_regression", params]]))) + @test_throws "needs ml_models_dir" resolve_aero_geometry( + yaml_in, joinpath(dir, "out"); verbose=false) + yaml_out = resolve_aero_geometry(yaml_in, joinpath(dir, "out"); + ml_models_dir=fixture_dir, verbose=false) + airfoil = YAML.load_file(yaml_out)["wing_airfoils"]["data"][1] + @test airfoil[2] == "polars" + section = Wing(yaml_out; n_panels=2).unrefined_sections[1] + @test section.aero_model == POLAR_VECTORS + cl, _, _ = masure_aero(load_masure_model(1e6, fixture_dir), params, [2.0]) + @test section.aero_data[1][4] ≈ deg2rad(2.0) + @test section.aero_data[2][4] ≈ cl[1] + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 39350918..9599109b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -70,6 +70,7 @@ function include_selected_tests() 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("airfoil_aero/test_masure_regression.jl") && include("airfoil_aero/test_masure_regression.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") # bin/release is a bash script, so only the unix runners can run it.