From 8fe1b8981e08b984899403e9794829e80268eb9d Mon Sep 17 00:00:00 2001 From: Anna Delbeke Date: Wed, 4 Mar 2026 10:45:22 +0100 Subject: [PATCH 01/10] Add wings3D problem template (copied from airfoil) --- engibench/problems/wings3D/README.md | 90 +++ engibench/problems/wings3D/__init__.py | 5 + .../problems/wings3D/dataset_slurm_airfoil.py | 167 ++++ .../wings3D/fake_pyoptsparse/__init__.py | 36 + .../fake_pyoptsparse/pyOpt_constraint.py | 6 + .../fake_pyoptsparse/pyOpt_objective.py | 6 + .../fake_pyoptsparse/pyOpt_optimization.py | 6 + .../fake_pyoptsparse/pyOpt_variable.py | 6 + engibench/problems/wings3D/pyopt_history.py | 740 ++++++++++++++++++ engibench/problems/wings3D/simulation_jobs.py | 56 ++ .../problems/wings3D/templates/__init__.py | 0 .../wings3D/templates/airfoil_analysis.py | 182 +++++ .../problems/wings3D/templates/airfoil_opt.py | 292 +++++++ .../wings3D/templates/cli_interface.py | 129 +++ .../problems/wings3D/templates/pre_process.py | 75 ++ engibench/problems/wings3D/utils.py | 429 ++++++++++ engibench/problems/wings3D/v0.py | 534 +++++++++++++ 17 files changed, 2759 insertions(+) create mode 100644 engibench/problems/wings3D/README.md create mode 100644 engibench/problems/wings3D/__init__.py create mode 100644 engibench/problems/wings3D/dataset_slurm_airfoil.py create mode 100644 engibench/problems/wings3D/fake_pyoptsparse/__init__.py create mode 100644 engibench/problems/wings3D/fake_pyoptsparse/pyOpt_constraint.py create mode 100644 engibench/problems/wings3D/fake_pyoptsparse/pyOpt_objective.py create mode 100644 engibench/problems/wings3D/fake_pyoptsparse/pyOpt_optimization.py create mode 100644 engibench/problems/wings3D/fake_pyoptsparse/pyOpt_variable.py create mode 100644 engibench/problems/wings3D/pyopt_history.py create mode 100644 engibench/problems/wings3D/simulation_jobs.py create mode 100644 engibench/problems/wings3D/templates/__init__.py create mode 100644 engibench/problems/wings3D/templates/airfoil_analysis.py create mode 100644 engibench/problems/wings3D/templates/airfoil_opt.py create mode 100644 engibench/problems/wings3D/templates/cli_interface.py create mode 100644 engibench/problems/wings3D/templates/pre_process.py create mode 100644 engibench/problems/wings3D/utils.py create mode 100644 engibench/problems/wings3D/v0.py diff --git a/engibench/problems/wings3D/README.md b/engibench/problems/wings3D/README.md new file mode 100644 index 00000000..353264a5 --- /dev/null +++ b/engibench/problems/wings3D/README.md @@ -0,0 +1,90 @@ +# Airfoil 2D + +**Lead**: Cashen Diniz @cashend + +Airfoil 2D is a benchmark problem that aims to optimize the shape of an airfoil to maximize the lift-to-drag ratio. +We rely on MACH-Aero for the simulations. + +## Side notes + +Here is the script I've used to upload the data to HF using the pickle files here: https://github.com/IDEALLab/OptimizingDiffusionSciTech2024/tree/main/data/optimized_data + +```python +from datasets import Dataset +from datasets import DatasetDict +import numpy as np +import pandas as pd + +opt_train_airfoils, opt_test_airfoils, opt_val_airfoils = pd.read_pickle("train_test_val_opt_airfoils.pkl") +init_train_airfoils, init_test_airfoils, init_val_airfoils = pd.read_pickle("train_test_val_init_airfoils.pkl") +train_params, test_params, val_params = pd.read_pickle("train_test_val_opt_params.pkl") + +# For each airfoil, we need one row containing the initial and optimized airfoil, as well as the parameters + +dataset_train = [] + +for o, i, p in zip(opt_train_airfoils, init_train_airfoils, train_params): + dataset_train.append( + { + "initial_design": {"coords": i, "angle_of_attack": np.asarray(p[4], dtype=np.float32)}, + "optimal_design": {"coords": o, "angle_of_attack": np.asarray(p[4], dtype=np.float32)}, + "mach": p[0], + "reynolds": p[1], + "cl_target": p[2], + "area_ratio_min": p[3], + "area_initial": p[5], + "cd": p[6], + "cl": p[7], + "cl_con_violation": p[8], + "area_ratio": p[9], + } + ) + +dataset_val = [] + +for o, i, p in zip(opt_test_airfoils, init_test_airfoils, test_params): + dataset_val.append( + { + "initial_design": {"coords": i, "angle_of_attack": np.asarray(p[4], dtype=np.float32)}, + "optimal_design": {"coords": o, "angle_of_attack": np.asarray(p[4], dtype=np.float32)}, + "mach": p[0], + "reynolds": p[1], + "cl_target": p[2], + "area_ratio_min": p[3], + "area_initial": p[5], + "cd": p[6], + "cl": p[7], + "cl_con_violation": p[8], + "area_ratio": p[9], + } + ) + +dataset_testt = [] + +for o, i, p in zip(opt_val_airfoils, init_val_airfoils, val_params): + dataset_testt.append( + { + "initial_design": {"coords": i, "angle_of_attack": np.asarray(p[4], dtype=np.float32)}, + "optimal_design": {"coords": o, "angle_of_attack": np.asarray(p[4], dtype=np.float32)}, + "mach": p[0], + "reynolds": p[1], + "cl_target": p[2], + "area_ratio_min": p[3], + "area_initial": p[5], + "cd": p[6], + "cl": p[7], + "cl_con_violation": p[8], + "area_ratio": p[9], + } + ) + + +# Create a huggingface dataset from the three splits above +train_spit = Dataset.from_list(dataset_train) +print(train_spit.shape) +val_spit = Dataset.from_list(dataset_val) +test_spit = Dataset.from_list(dataset_testt) +dataset_dict = DatasetDict({"train": train_spit, "val": val_spit, "test": test_spit}) +dataset_dict.push_to_hub("IDEALLab/airfoil_v0") + +``` diff --git a/engibench/problems/wings3D/__init__.py b/engibench/problems/wings3D/__init__.py new file mode 100644 index 00000000..1fae7ff8 --- /dev/null +++ b/engibench/problems/wings3D/__init__.py @@ -0,0 +1,5 @@ +"""Airfoil problem module.""" + +from engibench.problems.airfoil.v0 import Airfoil + +__all__ = ["Airfoil"] diff --git a/engibench/problems/wings3D/dataset_slurm_airfoil.py b/engibench/problems/wings3D/dataset_slurm_airfoil.py new file mode 100644 index 00000000..2eb05acc --- /dev/null +++ b/engibench/problems/wings3D/dataset_slurm_airfoil.py @@ -0,0 +1,167 @@ +"""Dataset Generation for Airfoil Problem via SLURM. + +This script generates a dataset for the Airfoil problem using the SLURM API +""" + +from argparse import ArgumentParser + +from datasets import load_dataset +import numpy as np +from scipy.stats import qmc + +from engibench.problems.airfoil.simulation_jobs import simulate_slurm +from engibench.utils import slurm + + +def calculate_runtime(group_size, minutes_per_sim=5): + """Calculate runtime based on group size and (rough) estimate of minutes per simulation.""" + total_minutes = group_size * minutes_per_sim + hours = total_minutes // 60 + minutes = total_minutes % 60 + return f"{hours:02d}:{minutes:02d}:00" + + +if __name__ == "__main__": + """Dataset Generation, Simulation, and Rendering for Airfoil Problem via SLURM. + + This script generates a dataset for the Airfoil problem using the SLURM API, though it could + be generalized to other problems as well. It includes functions for simulation of designs. + + Command Line Arguments: + -n_designs, --num_designs: How many airfoil designs should we use? + -n_flows, --num_flow_conditions: How many flow conditions should we use per design? + -n_aoas, --num_angles_of_attack: How many angles of attack should we use per design & flow condition pairing? + -group_size, --group_size: How many simulations should we group together on a single cpu? + -n_slurm_array, --num_slurm_array: How many slurm jobs to spawn and submit via slurm arrays? Note this may be limited by the HPC system. + """ + # Fetch command line arguments for render and simulate to know whether to run those functions + parser = ArgumentParser() + parser.add_argument( + "-n_designs", + "--num_designs", + type=int, + default=10, + help="How many airfoil designs should we use?", + ) + parser.add_argument( + "-n_flows", + "--num_flow_conditions", + type=int, + default=1, + help="How many flow conditions (Mach Number and Reynolds Number) should we sample for each design?", + ) + parser.add_argument( + "-n_aoas", + "--num_angles_of_attack", + type=int, + default=1, + help="How many angles of attack should we sample for each design?", + ) + parser.add_argument( + "-group_size", + "--group_size", + type=int, + default=2, + help="How many simulations do you wish to batch within each individual slurm job?", + ) + parser.add_argument( + "-n_slurm_array", + "--num_slurm_array", + type=int, + default=1000, + help="What is the maximum size of the Slurm array (Will vary from HPC system to HPC system)?", + ) + args = parser.parse_args() + + n_designs = args.num_designs + n_flows = args.num_flow_conditions + n_aoas = args.num_angles_of_attack + group_size = args.group_size + n_slurm_array = args.num_slurm_array + + # ============== Problem-specific elements =================== + # The following elements are specific to the problem and should be modified accordingly + + # Define flow parameter and angle of attack ranges + Ma_min, Ma_max = 0.5, 0.9 # Mach number range + Re_min, Re_max = 1.0e6, 2.0e7 # Reynolds number range + aoa_min, aoa_max = 0.0, 20.0 # Angle of attack range + + # Load airfoil designs from HF Database + ds = load_dataset("IDEALLab/airfoil_v0") + designs = ( + ds["train"]["initial_design"] + + ds["train"]["optimal_design"] + + ds["val"]["initial_design"] + + ds["val"]["optimal_design"] + + ds["test"]["initial_design"] + + ds["test"]["optimal_design"] + ) + + # Use specified number of designs + designs = designs[:n_designs] + + # Generate LHS samples + rng = np.random.default_rng(seed=42) # Optional seed for reproducibility + sampler = qmc.LatinHypercube(d=2, seed=rng) + samples = sampler.random(n=n_designs * n_flows) # n samples needed + + # Scale to your flow domain + bounds = np.array([[Ma_min, Ma_max], [Re_min, Re_max]]) + scaled_samples = qmc.scale(samples, bounds[:, 0], bounds[:, 1]) + mach_values = scaled_samples[:, 0] + reynolds_values = scaled_samples[:, 1] + + # Generate all simulation configurations + config_id = 0 + simulate_configs_designs = [] + for i, design in enumerate(designs): + for j in range(n_flows): + ma = mach_values[i * n_flows + j] + re = reynolds_values[i * n_flows + j] + for alpha in rng.uniform(low=aoa_min, high=aoa_max, size=n_aoas): + problem_configuration = {"mach": ma, "reynolds": re, "alpha": alpha} + config = {"problem_configuration": problem_configuration, "configuration_id": config_id} + config["design"] = design["coords"] + simulate_configs_designs.append(config) + config_id += 1 + + print(f"Generated {len(simulate_configs_designs)} configurations for simulation.") + + # Calculate total number of simulation jobs and number of sbatch maps needed + n_simulations = len(simulate_configs_designs) + n_sbatch_maps = np.ceil(n_simulations / (group_size * n_slurm_array)) + + slurm_config = slurm.SlurmConfig( + name="Airfoil_dataset_generation", + runtime=calculate_runtime(group_size, minutes_per_sim=5), + ntasks=1, + cpus_per_task=1, + log_dir="./sim_logs/", + ) + print(calculate_runtime(group_size, minutes_per_sim=5)) + + submitted_jobs = [] + for ibatch in range(int(n_sbatch_maps)): + sim_batch_configs = simulate_configs_designs[ + ibatch * group_size * n_slurm_array : (ibatch + 1) * group_size * n_slurm_array + ] + print(len(sim_batch_configs)) + print(f"Submitting batch {ibatch + 1}/{int(n_sbatch_maps)}") + + job_array = slurm.sbatch_map( + f=simulate_slurm, + args=sim_batch_configs, + slurm_args=slurm_config, + group_size=group_size, # Number of jobs to batch in sequence to reduce job array size + work_dir="scratch", + ) + + # Save the job array reference + submitted_jobs.append(job_array) + + # Wait for this job to complete by calling save() + # This will submit a dependent job that waits for the array to finish + print(f"Waiting for batch {ibatch + 1} to complete...") + job_array.save(f"results_{ibatch}.pkl", slurm_args=slurm_config) + print(f"Batch {ibatch + 1} completed!") diff --git a/engibench/problems/wings3D/fake_pyoptsparse/__init__.py b/engibench/problems/wings3D/fake_pyoptsparse/__init__.py new file mode 100644 index 00000000..e05bd8d8 --- /dev/null +++ b/engibench/problems/wings3D/fake_pyoptsparse/__init__.py @@ -0,0 +1,36 @@ +"""Drop-in module for pyoptsparse to unpickle ahistory when pyoptsparse is not installed.""" + +from types import ModuleType + + +class FakePyOptSparseObject: + """Drop-in for objects needed to unpickle a pyoptsparse history when pyoptsparse is not installed.""" + + def __init__(self, *args, **kwargs) -> None: + self.args = args + self.kwargs = kwargs + + def _mapContoOpt_Dict(self, d): # noqa: N802 + return d + + def _mapXtoOpt_Dict(self, d): # noqa: N802 + return d + + def _mapObjtoOpt_Dict(self, d): # noqa: N802 + return d + + +class Optimization(FakePyOptSparseObject): + """Drop-in.""" + + +class Variable(FakePyOptSparseObject): + """Drop-in.""" + + +class Constraint(FakePyOptSparseObject): + """Drop-in.""" + + +class Objective(FakePyOptSparseObject): + """Drop-in.""" diff --git a/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_constraint.py b/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_constraint.py new file mode 100644 index 00000000..5e006bf9 --- /dev/null +++ b/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_constraint.py @@ -0,0 +1,6 @@ +# noqa: N999 +"""Drop-in.""" + +from engibench.problems.airfoil.fake_pyoptsparse import FakePyOptSparseObject as Constraint + +__all__ = ["Constraint"] diff --git a/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_objective.py b/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_objective.py new file mode 100644 index 00000000..8df45980 --- /dev/null +++ b/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_objective.py @@ -0,0 +1,6 @@ +# noqa: N999 +"""Drop-in.""" + +from engibench.problems.airfoil.fake_pyoptsparse import FakePyOptSparseObject as Objective + +__all__ = ["Objective"] diff --git a/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_optimization.py b/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_optimization.py new file mode 100644 index 00000000..b8eadd60 --- /dev/null +++ b/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_optimization.py @@ -0,0 +1,6 @@ +# noqa: N999 +"""Drop-in.""" + +from engibench.problems.airfoil.fake_pyoptsparse import FakePyOptSparseObject as Optimization + +__all__ = ["Optimization"] diff --git a/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_variable.py b/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_variable.py new file mode 100644 index 00000000..ec00fc6f --- /dev/null +++ b/engibench/problems/wings3D/fake_pyoptsparse/pyOpt_variable.py @@ -0,0 +1,6 @@ +# noqa: N999 +"""Drop-in.""" + +from engibench.problems.airfoil.fake_pyoptsparse import FakePyOptSparseObject as Variable + +__all__ = ["Variable"] diff --git a/engibench/problems/wings3D/pyopt_history.py b/engibench/problems/wings3D/pyopt_history.py new file mode 100644 index 00000000..ef0fb76b --- /dev/null +++ b/engibench/problems/wings3D/pyopt_history.py @@ -0,0 +1,740 @@ +# ruff: noqa: D101,D205,N803,N802,N806,FBT002,D416,SIM118,RET503,RET502,RET505,D212,ISC003,SLF001,C901,PLR0913,PLR0912,PLR0915,S110,BLE001,SIM102,PLR2004 +"""This file is copy pasted from pyOptSparse: https://github.com/mdolab/pyoptsparse/blob/main/pyoptsparse/pyOpt_history.py. + +It is used to read the history of an optimization run. +""" + +# Standard Python modules +from collections import OrderedDict +import copy +import os + +# External modules +import numpy as np +from sqlitedict import SqliteDict + +EPS = np.finfo(np.float64).eps + + +class pyOptSparseWarning: # noqa: N801 + """Format a warning message.""" + + def __init__(self, message): + msg = "\n+" + "-" * 78 + "+" + "\n" + "| pyOptSparse Warning: " + i = 23 + for word in message.split(): + if len(word) + i + 1 > 78: # Finish line and start new one + msg += " " * (79 - i) + "|\n| " + word + " " + i = 1 + len(word) + 1 + else: + msg += word + " " + i += len(word) + 1 + msg += " " * (78 - i) + "|\n" + "+" + "-" * 78 + "+" + "\n" + print(msg) + + +class History: + def __init__(self, fileName, optProb=None, temp=False, flag="r"): + """This class is essentially a thin wrapper around a SqliteDict dictionary to facilitate + operations with pyOptSparse. + + Parameters + ---------- + fileName : str + File name for history file + + optProb : pyOpt_Optimization + The optimization object + + temp : bool + Flag to signify that the file should be deleted after it is + closed + + flag : str + String specifying the mode. Similar to what was used in shelve. + ``n`` for a new database and ``r`` to read an existing one. + """ + self.flag = flag + if self.flag == "n": + # If writing, we expliclty remove the file to + # prevent old keys from "polluting" the new histrory + if os.path.exists(fileName): + os.remove(fileName) + self.db = SqliteDict(fileName) + self.optProb = optProb + elif self.flag == "r": + if os.path.exists(fileName): + # we cast the db to OrderedDict so we do not have to + # manually close the underlying db at the end + self.db = OrderedDict(SqliteDict(fileName)) + else: + raise FileNotFoundError(f"The requested history file {fileName} to open in read-only mode does not exist.") + self._processDB() + else: + raise ValueError("The flag argument to History must be 'r' or 'n'.") + self.temp = temp + self.fileName = fileName + + def close(self): + """Close the underlying database. + This should only be used in write mode. In read mode, we close the db + during initialization. + """ + if self.flag == "n": + self.db.close() + if self.temp: + os.remove(self.fileName) + + def write(self, callCounter, data): + """This is the main to write data. Basically, we just pass in + the callCounter, the integer forming the key, and a dictionary + which will be written to the key. + + Parameters + ---------- + callCounter : int + the callCounter to write to + data : dict + the dictionary corresponding to the callCounter + """ + # String key to database on disk + key = str(callCounter) + # if the point exists, we merely update with new data + if self.pointExists(callCounter): + oldData = self.read(callCounter) + oldData.update(data) + self.db[key] = oldData + else: + self.db[key] = data + self.db["last"] = key + self.db.sync() + self.keys = list(self.db.keys()) + + def writeData(self, key, data): + """Write arbitrary `key:data` value to db. + + Parameters + ---------- + key : str + The key to be added to the history file + data + The data corresponding to the key. It can be anything as long as it is serializable + in `sqlitedict`. + """ + self.db[key] = data + self.db.commit() + self.keys = list(self.db.keys()) + + def pointExists(self, callCounter): + """Determine if callCounter is in the database. + + Parameters + ---------- + callCounter : int or str of int + + Returns + ------- + bool + True if the callCounter exists in the history file. + False otherwise. + """ + if isinstance(callCounter, int): + callCounter = str(callCounter) + return callCounter in self.keys + + def read(self, key): + """Read data for an arbitrary key. Returns None if key is not found. + If passing in a callCounter, it should be verified by calling pointExists() first. + + Parameters + ---------- + key : str or int + generic key[str] or callCounter[int] + + Returns + ------- + dict + The value corresponding to `key` is returned. + If the key is not found, `None` is returned instead. + """ + if isinstance(key, int): + key = str(key) + try: + return self.db[key] + except KeyError: + return None + + def _searchCallCounter(self, x): + """Searches through existing callCounters, and finds the one corresponding + to an evaluation at the design vector `x`. + returns `None` if the point did not match previous evaluations. + + Parameters + ---------- + x : ndarray + The unscaled DV as a single array. + + Returns + ------- + int + The callCounter corresponding to the DV `x`. + `None` is returned if no match was found. + + Notes + ----- + The tolerance used for this is the value `numpy.finfo(numpy.float64).eps`. + """ + last = int(self.db["last"]) + callCounter = None + for i in range(last, 0, -1): + key = str(i) + xuser = self.optProb.processXtoVec(self.db[key]["xuser"]) + if np.isclose(xuser, x, atol=EPS, rtol=EPS).all() and "funcs" in self.db[key].keys(): + callCounter = i + break + return callCounter + + def _processDB(self): + """Pre-processes the DB file and store various values into class attributes. + These will be used later when calling self.getXX functions. + """ + # Load any keys it happens to have: + self.keys = list(self.db.keys()) + # load info + self.DVInfo = self.read("varInfo") + self.conInfo = self.read("conInfo") + self.objInfo = self.read("objInfo") + # metadata + self.metadata = self.read("metadata") + self.optProb = self.read("optProb") + # load names + self.DVNames = set(self.getDVNames()) + self.conNames = set(self.getConNames()) + self.objNames = set(self.getObjNames()) + + # extract list of callCounters from self.keys + # this just checks if each key contains only digits, then cast into int + self.callCounters = sorted([x for x in self.keys if x.isdigit()], key=float) + + # extract all information stored in the call counters + self.iterKeys = set() + self.extraFuncsNames = set() + for i in self.callCounters: + val = self.read(i) + self.iterKeys.update(val.keys()) + if "funcs" in val.keys(): + self.extraFuncsNames.update(val["funcs"].keys()) + # remove objective and constraint keys + self.extraFuncsNames = self.extraFuncsNames.difference(self.conNames).difference(self.objNames) + + if self.metadata["version"] != "2.13.2": + pyOptSparseWarning( + "The version of pyoptsparse used to generate the history file does not match the one being run right now. There may be compatibility issues." + ) + + def getIterKeys(self): + """Returns the keys available at each optimization iteration. + This function is useful for inspecting the history file, to determine + what information is saved at each iteration. + + Returns + ------- + list of str + A list containing the names of keys stored at each optimization iteration. + """ + return copy.deepcopy(list(self.iterKeys)) + + def getDVNames(self): + """Returns the names of the DVs. + + Returns + ------- + list of str + A list containing the names of DVs. + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + return copy.deepcopy(list(self.DVInfo.keys())) + + def getConNames(self): + """Returns the names of constraints. + + Returns + ------- + list of str + A list containing the names of constraints. + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + # we remove linear constraints + conNames = [con for con in self.conInfo.keys() if not self.optProb.constraints[con].linear] + return copy.deepcopy(conNames) + + def getObjNames(self): + """Returns the names of the objectives. + + Returns + ------- + list of str + A list containing the names of objectives. + + Notes + ----- + Recall that for the sake of generality, pyOptSparse allows for multiple objectives to be + added. This feature is not used currently, but does make `ObjNames` follow the same structure + as `ConNames` and `DVNames`. + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + return copy.deepcopy(list(self.objInfo.keys())) + + def getExtraFuncsNames(self): + """Returns extra funcs names. + These are extra key: value pairs stored in the ``funcs`` dictionary for each iteration, which are not used by the optimizer. + + Returns + ------- + list of str + A list containing the names of extra funcs keys. + + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + return copy.deepcopy(list(self.extraFuncsNames)) + + def getObjInfo(self, key=None): + """Returns the `ObjInfo`, for all keys by default. A `key` parameter can also + be supplied, to retrieve `ObjInfo` corresponding to specific keys. + + + Parameters + ---------- + key : str or list of str, optional + Specifies for which obj to extract `ObjInfo`. + + Returns + ------- + dict + A dictionary containing ObjInfo. For a single key, the return is one level deeper. + + Notes + ----- + Recall that for the sake of generality, pyOptSparse allows for multiple objectives to be + added. This feature is not used currently, but does make `ObjInfo` follow the same structure + as `ConInfo` and `DVInfo`. + Because of this, it is recommended that this function be accessed using the optional `key` + argument. + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + if key is not None: + if isinstance(key, str): + return copy.deepcopy(self.objInfo[key]) + elif isinstance(key, list): + d = OrderedDict() + for k in key: + d[k] = self.objInfo[k] + return d + else: + return copy.deepcopy(self.objInfo) + + def getConInfo(self, key=None): + """Returns the `ConInfo`, for all keys by default. A `key` parameter can also + be supplied, to retrieve `ConInfo` corresponding to specific constraints. + + Parameters + ---------- + key : str or list of str, optional + Specifies for which constraint to extract `ConInfo`. + + Returns + ------- + dict + A dictionary containing ConInfo. For a single key, the return is one level deeper. + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + if key is not None: + if isinstance(key, str): + return copy.deepcopy(self.conInfo[key]) + elif isinstance(key, list): + d = OrderedDict() + for k in key: + d[k] = self.conInfo[k] + return d + else: + return copy.deepcopy(self.conInfo) + + def getDVInfo(self, key=None): + """Returns the `DVInfo`, for all keys by default. A `key` parameter can also + be supplied, to retrieve `DVInfo` corresponding to specific DVs. + + Parameters + ---------- + key : str or list of str, optional + Specifies for which DV to extract `DVInfo`. + + Returns + ------- + dict + A dictionary containing DVInfo. For a single key, the return is one level deeper. + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + if key is not None: + if isinstance(key, str): + return copy.deepcopy(self.DVInfo[key]) + elif isinstance(key, list): + d = OrderedDict() + for k in key: + d[k] = self.DVInfo[k] + return d + else: + return copy.deepcopy(self.DVInfo) + + def getMetadata(self): + """ + Returns a copy of the metadata stored in the history file. + + Returns + ------- + dict + A dictionary containing the metadata. + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + return copy.deepcopy(self.metadata) + + def getOptProb(self): + """ + Returns a copy of the optProb associated with the optimization. + + Returns + ------- + optProb : pyOpt_optimization object + The optProb associated with the optimization. This is taken from the history file, + and therefore has the ``comm`` and ``objFun`` fields removed. + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + return copy.deepcopy(self.optProb) + + def _processIterDict(self, d, scale=False): + """ + This function scales the value, where the factor is extracted from the + `Info` dictionaries, according to "name". + + Parameters + ---------- + d : dictionary + The iteration dictionary, i.e. hist['0'] + This must be a function evaluation callCounter, and not a gradient callCounter. + scale : bool + Whether the returned values should be scaled. + + Returns + ------- + conDict : dict + A dictionary containing constraint values + objDict : dict + A dictionary containing objective values + DVDict : dict + A dictionary containing DV values + + These are all "flat" dictionaries, with simple key:value pairs. + """ + conDict = {} + objDict = {} + # these require funcs which may not always be there + if "funcs" in d.keys(): + for con in list(self.optProb.constraints.keys()): + # linear constraints are not stored in funcs + if not self.optProb.constraints[con].linear: + conDict[con] = d["funcs"][con] + else: + # the linear constraints are removed from optProb so that scaling works + # without needing the linear constraints to be present + self.optProb.constraints.pop(con) + + for obj in self.objNames: + objDict[obj] = d["funcs"][obj] + + # the DVs will always be there + DVDict = {} + for DV in self.DVNames: + DVDict[DV] = d["xuser"][DV] + if scale: + conDict = self.optProb._mapContoOpt_Dict(conDict) + objDict = self.optProb._mapObjtoOpt_Dict(objDict) + DVDict = self.optProb._mapXtoOpt_Dict(DVDict) + return conDict, objDict, DVDict + + def getCallCounters(self): + """ + Returns a list of all call counters stored in the history file. + + Returns + ------- + list + a list of strings, each entry being a call counter. + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + return copy.deepcopy(self.callCounters) + + def getValues(self, names=None, callCounters=None, major=True, scale=False, stack=False, allowSens=False): + """ + Parses an existing history file and returns a data dictionary used to post-process optimization results, containing the requested optimization iteration history. + + Parameters + ---------- + names : list or str + the values of interest, can be the name of any DV, objective or constraint, + or a list of them. If None, all values are returned. This includes the DVs, + funcs, and any values stored by the optimizer. + + callCounters : list of ints, can also contain 'last' + a list of callCounters to extract information from. + If the callCounter is invalid, i.e. outside the range or is a funcsSens evaluation, then it is skipped. + 'last' represents the last major iteration. + If None, values from all callCounters are returned. + + major : bool + flag to specify whether to include only major iterations. + + scale : bool + flag to specify whether to apply scaling for the values. True means + to return values that are scaled the same way as the actual optimization. + + stack : bool + flag to specify whether the DV should be stacked into a single numpy array with + the key `xuser`, or retain their separate DVGroups. + + allowSens: bool + flag to specify whether gradient evaluation iterations are allowed. + If true, it is up to the user to ensure that the callCounters specified contain the information requested. + + Returns + ------- + dict + a dictionary containing the information requested. The keys of the dictionary + correspond to the `names` requested. Each value is a numpy array with the first + dimension equal to the number of callCounters requested. + + Notes + ----- + Regardless of the major flag, failed function evaluations are not returned. + + Examples + -------- + First we can request DV history over all major iterations: + + >>> hist.getValues(names="xvars", major=True) + {'xvars': array([[-2.00000000e+00, 1.00000000e+00], + [-1.00000000e+00, 9.00000000e-01], + [-5.00305827e-17, 4.21052632e-01], + [ 1.73666171e-06, 4.21049838e-01], + [ 9.08477459e-06, 4.21045261e-01], + [ 5.00000000e-01, 2.84786405e-01], + [ 5.00000000e-01, 5.57279939e-01], + [ 5.00000000e-01, 2.00000000e+00]])} + + Next we can look at DV and optimality for the first and last iteration only: + + >>> hist.getValues(names=["xvars", "optimality"], callCounters=[0, "last"]) + {'optimality': array([1.27895528, 0. ]), + 'xvars': array([[-2. , 1. ], + [ 0.5, 2. ]])} + """ + # only do this if we open the file with 'r' flag + if self.flag != "r": + return + + allNames = ( + self.DVNames.union(self.conNames) + .union(self.objNames) + .union(self.iterKeys) + .union(self.extraFuncsNames) + .difference({"funcs", "funcsSens", "xuser"}) + ) + # cast string input into a single list + if isinstance(names, str): + names = {names} + elif names is None: + names = allNames + else: + names = set(names) + if stack: + allNames.add("xuser") + # error if names isn't either a DV, con or obj + if not names.issubset(allNames): + raise KeyError( + "The names provided are not one of DVNames, conNames or objNames.\n" + + f"The names must be a subset of {allNames}" + ) + DVsAsFuncs = self.DVNames.intersection(self.conNames) + if len(DVsAsFuncs) > 0: + ambiguousNames = names.intersection(DVsAsFuncs) + if len(ambiguousNames) > 0: + pyOptSparseWarning( + f"The names provided {ambiguousNames} is ambiguous, since it is both a DV as well as an objective/constraint. " + + "It is being assumed to be a DV. If it was set up via addDVsAsFunctions, then there's nothing to worry. " + + "Otherwise, consider renaming the variable or manually editing the history file." + ) + + if len(names.intersection(self.iterKeys)) > 0: + if not major: + pyOptSparseWarning( + "The major flag has been set to True, since some names specified only exist on major iterations." + ) + major = True + + if stack: + DVinNames = names.intersection(self.DVNames) + for DV in DVinNames: + names.remove(DV) + names.add("xuser") + pyOptSparseWarning( + "The stack flag was set to True. Therefore all DV names have been removed, and replaced with a single key 'xuser'." + ) + + # set up dictionary to return + data = {} + # pre-allocate list for each input + for name in names: + data[name] = [] + + # this flag is used for error printing only + user_specified_callCounter = False + if callCounters is not None: + user_specified_callCounter = True + if isinstance(callCounters, str): + callCounters = [callCounters] + else: + callCounters = self.callCounters + + # parse the 'last' callCounter + if "last" in callCounters: + callCounters.append(self.read("last")) + callCounters.remove("last") + + self._previousIterCounter = -1 + # loop over call counters, check if each counter is valid, and parse + for i in callCounters: + val = self._readValidCallCounter(i, user_specified_callCounter, allowSens, major) + if val is not None: # if i is valid + conDict, objDict, DVDict = self._processIterDict(val, scale=scale) + for name in names: + if name == "xuser": + data[name].append(self.optProb.processXtoVec(DVDict)) + elif name in self.DVNames: + data[name].append(DVDict[name]) + elif name in self.conNames: + data[name].append(conDict[name]) + elif name in self.objNames: + data[name].append(objDict[name]) + elif name in self.extraFuncsNames: + data[name].append(val["funcs"][name]) + else: # must be opt + data[name].append(val[name]) + + # reshape lists into numpy arrays + for name in names: + # we just stack along axis 0 + if len(data[name]) > 0: + data[name] = np.stack(data[name], axis=0) + else: + data[name] = np.array(data[name]) + # we cast 1D arrays to 2D, for scalar values + if data[name].ndim == 1: + data[name] = np.expand_dims(data[name], 1) + + # Raise warning for IPOPT's duplicated history + if self.db["metadata"]["optimizer"] == "IPOPT" and "iter" not in self.db["0"].keys(): + pyOptSparseWarning( + "The optimization history of IPOPT has duplicated entries at every iteration. " + + "Fix the history manually, or re-run the optimization with a current version of pyOptSparse to generate a correct history file. " + ) + return data + + def _readValidCallCounter(self, i, user_specified_callCounter, allowSens, major): + """ + Checks whether a call counter is valid and read the data. The call counter is valid when it is + 1) inside the range of the history data, + 2) a function evaluation (i.e. not a sensitivity evaluation, except when `allowSens = True`), + 3) not a duplicated entry, + 4) not a failed function evaluation, + 5) a major iteration (only when `major = True`). + + Parameters + ---------- + i : int + call counter. + + user_specified_callCounter : bool + flag to specify whether the call counter `i` is requested by a user or not. + + allowSens: bool + flag to specify whether gradient evaluation iterations are allowed. + + major : bool + flag to specify whether to include only major iterations. + + Returns + ------- + val : dict or None + information corresponding to the call counter `i`. + If the call counter is not valid, `None` is returned instead. + """ + if not self.pointExists(i): + if user_specified_callCounter: + # user specified a non-existent call counter + pyOptSparseWarning(f"callCounter {i} was not found and is skipped!") + return None + else: + val = self.read(i) + + # check if the callCounter is of a function call + if not ("funcs" in val.keys() or allowSens): + if user_specified_callCounter: + # user unintentionally specified a call counter for sensitivity + pyOptSparseWarning( + f"callCounter {i} did not contain a function evaluation and is skipped! " + + "Was it a gradient evaluation step?" + ) + return None + else: + # exclude the duplicated history (only when we have "iter" recorded) + if "iter" in val.keys(): + duplicate_flag = val["iter"] == self._previousIterCounter + self._previousIterCounter = val["iter"] # update iterCounter for next i + if duplicate_flag and not user_specified_callCounter: + # this is a duplicate + return None + # end if "iter" in val.keys() + + # check major/minor iteration, and if the call failed + if ((major and val["isMajor"]) or not major) and not val["fail"]: + return val + else: + return None + # end if - ("funcs" in val.keys() + # end if - pointExists + + def __del__(self): + try: + self.db.close() + if self.temp: + os.remove(self.fileName) + except Exception: + pass diff --git a/engibench/problems/wings3D/simulation_jobs.py b/engibench/problems/wings3D/simulation_jobs.py new file mode 100644 index 00000000..265b5410 --- /dev/null +++ b/engibench/problems/wings3D/simulation_jobs.py @@ -0,0 +1,56 @@ +"""Dataset Generator for the Airfoil problem using the SLURM API.""" + +import time + +import numpy as np + +from engibench.problems.airfoil.v0 import Airfoil + + +def simulate_slurm(problem_configuration: dict, configuration_id: int, design: list) -> dict: + """Takes in the given configuration and designs, then runs the simulation analysis. + + Any arguments should be things that you want to change across the different jobs, and anything + that is the same/static across the runs should just be defined inside this function. + + Args: + problem_configuration (dict): The specific configuration used to setup the problem being passed. + For the airfoil problem this includes Mach number, Reynolds number, and angle of attack. + configuration_id (int): A unique identifier for the job for later debugging or tracking. + design (list): list of lists defining x and y coordinates of airfoil geometry. + + Returns: + "performance_dict": Dictionary of aerodynamic performance (lift & drag). + "simulate_time": The time taken to run this simulation job. Useful for aggregating + the time taken for dataset generation. + "problem_configuration": Problem configuration parameters + "configuration_id": Identifier for specific simulation configurations + """ + # Instantiate problem + problem = Airfoil() + + # Set simulation ID + sim_id = configuration_id + 1 + + # Create unique simulation directory + problem.reset(seed=sim_id, cleanup=False) + + # Create simulation design (coordinates + angle of attack) + my_design = {"coords": np.array(design), "angle_of_attack": problem_configuration["alpha"]} + + print("Starting `simulate` via SLURM...") + start_time = time.time() + + performance = problem.simulate(my_design, mpicores=1, config=problem_configuration) + performance_dict = {"drag": performance[0], "lift": performance[1]} + print("Finished `simulate` via SLURM.") + end_time = time.time() + elapsed_time = end_time - start_time + print(f"Elapsed time for `simulate`: {elapsed_time:.2f} seconds") + + return { + "performance_dict": performance_dict, + "simulate_time": elapsed_time, + "problem_configuration": problem_configuration, + "configuration_id": configuration_id, + } diff --git a/engibench/problems/wings3D/templates/__init__.py b/engibench/problems/wings3D/templates/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/engibench/problems/wings3D/templates/airfoil_analysis.py b/engibench/problems/wings3D/templates/airfoil_analysis.py new file mode 100644 index 00000000..e73b1ebd --- /dev/null +++ b/engibench/problems/wings3D/templates/airfoil_analysis.py @@ -0,0 +1,182 @@ +"""This file is based on the MACHAero tutorials. + +https://github.com/mdolab/MACH-Aero/blob/main/tutorial/ +""" + +import itertools +import json +import os +import sys +from typing import Any + +from adflow import ADFLOW +from baseclasses import AeroProblem +from cli_interface import AnalysisParameters # type:ignore[import-not-found] +from cli_interface import Task +from mpi4py import MPI +import numpy as np + + +def main() -> None: # noqa: C901, PLR0915 + """Entry point of the script.""" + args = AnalysisParameters.from_dict(json.loads(sys.argv[1])) + mesh_fname = args.mesh_fname + output_dir = args.output_dir + task = args.task + + # mach number + mach = args.mach + # Reynolds number + reynolds = args.reynolds + # altitude + altitude = args.altitude + # temperature + temperature = args.temperature + # Whether to use altitude + use_altitude = args.use_altitude + # Reynold's Length + reynolds_length = 1.0 + + comm = MPI.COMM_WORLD + print(f"Processor {comm.rank} of {comm.size} is running") + if not os.path.exists(output_dir) and comm.rank == 0: + os.mkdir(output_dir) + + # rst ADflow options + aero_options = { + # I/O Parameters + "gridFile": mesh_fname, + "outputDirectory": output_dir, + "monitorvariables": ["cl", "cd", "resrho", "resrhoe"], + "writeTecplotSurfaceSolution": True, + # Physics Parameters + "equationType": "RANS", + "smoother": "DADI", + "rkreset": True, + "nrkreset": 10, + # Solver Parameters + "MGCycle": "sg", + # ANK Solver Parameters + "useANKSolver": True, + "ankswitchtol": 1e-1, + "liftIndex": 2, + "nsubiterturb": 10, + # NK Solver Parameters + "useNKSolver": True, + "NKSwitchTol": 1e-4, + # Termination Criteria + "L2Convergence": 1e-9, + "L2ConvergenceCoarse": 1e-4, + "nCycles": 5000, + } + print("ADflow options:") + # rst Start ADflow + # Create solver + cfd_solver = ADFLOW(options=aero_options) + + # Add features + span = 1.0 + pos = np.array([0.5]) * span + cfd_solver.addSlices("z", pos, sliceType="absolute") + + # rst Create AeroProblem + alpha = args.alpha + + if use_altitude: + ap = AeroProblem( + name="fc", + alpha=alpha, + mach=mach, + altitude=altitude, + areaRef=1.0, + chordRef=1.0, + evalFuncs=["cl", "cd"], + ) + else: + ap = AeroProblem( + name="fc", + alpha=alpha, + mach=mach, + T=temperature, + reynolds=reynolds, + reynoldsLength=reynolds_length, + areaRef=1.0, + chordRef=1.0, + evalFuncs=["cl", "cd"], + ) + + # rst Run ADflow + if task == Task.ANALYSIS: + print("Running analysis") + # Solve + cfd_solver(ap) + # rst Evaluate and print + funcs: dict[str, Any] = {} + cfd_solver.evalFunctions(ap, funcs) + # Print the evaluated functions + if comm.rank == 0: + cl = funcs[f"{ap.name}_cl"] + cd = funcs[f"{ap.name}_cd"] + # Save the lift and drag coefficients to a file + outputs = np.array([mach, reynolds, alpha, cl, cd]) + np.save(os.path.join(output_dir, "outputs.npy"), outputs) + + # rst Create polar arrays + elif task == Task.POLAR: + print("Running polar") + # Create an array of alpha values. + # In this case we create 5 random alpha values between 0 and 10 + alphas = np.linspace(0, 20, 50) + # Sort the alpha values + alphas.sort() + + # Create storage for the evaluated lift and drag coefficients + cl_list = [] + cd_list = [] + reslist = [] + # rst Start loop + # Loop over the alpha values and evaluate the polar + for alpha in alphas: + # rst update AP + # Update the name in the AeroProblem. This allows us to modify the + # output file names with the current alpha. + ap.name = f"fc_{alpha:4.2f}" + + # Update the alpha in aero problem and print it to the screen. + ap.alpha = alpha + if comm.rank == 0: + print(f"current alpha: {ap.alpha}") + + # rst Run ADflow polar + # Solve the flow + cfd_solver(ap) + + # Evaluate functions + funcs = {} + cfd_solver.evalFunctions(ap, funcs) + + # Store the function values in the output list + cl_list.append(funcs[f"{ap.name}_cl"]) + cd_list.append(funcs[f"{ap.name}_cd"]) + reslist.append(cfd_solver.getFreeStreamResidual(ap)) + if comm.rank == 0: + print(f"CL: {cl_list[-1]}, CD: {cd_list[-1]}") + + # rst Print polar + # Print the evaluated functions in a table + if comm.rank == 0: + outputs = np.array( + list( + zip(itertools.repeat(mach), itertools.repeat(reynolds), alphas, cl_list, cd_list, reslist, strict=False) + ) + ) + for *_, alpha_v, cl, cd, _res in outputs: + print(f"{alpha_v:6.1f} {cl:8.4f} {cd:8.4f}") + # Save the lift and drag coefficients to a file + np.save(os.path.join(output_dir, "M_Re_alpha_CL_CD_res.npy"), outputs) + + MPI.COMM_WORLD.Barrier() + + +if __name__ == "__main__": + main() diff --git a/engibench/problems/wings3D/templates/airfoil_opt.py b/engibench/problems/wings3D/templates/airfoil_opt.py new file mode 100644 index 00000000..6a0b9c75 --- /dev/null +++ b/engibench/problems/wings3D/templates/airfoil_opt.py @@ -0,0 +1,292 @@ +"""This file is largely based on the MACHAero tutorials. + +https://github.com/mdolab/MACH-Aero/blob/main/tutorial/ +""" + +# ====================================================================== +# Import modules +# ====================================================================== +import json +import os +import sys + +from adflow import ADFLOW +from baseclasses import AeroProblem +from cli_interface import Algorithm # type:ignore[import-not-found] +from cli_interface import OptimizeParameters +from idwarp import USMesh +from mpi4py import MPI +from multipoint import multiPointSparse +import numpy as np +from pygeo import DVConstraints +from pygeo import DVGeometry +from pyoptsparse import OPT +from pyoptsparse import Optimization + + +def main() -> None: # noqa: C901, PLR0915 + """Entry point of the script.""" + args = OptimizeParameters.from_dict(json.loads(sys.argv[1])) + # ====================================================================== + # Specify parameters for optimization + # ====================================================================== + # cL constraint + mycl = args.cl_target + # angle of attack + alpha = args.alpha + # mach number + mach = args.mach + # Reynolds number + reynolds = args.reynolds + # cruising altitude + altitude = args.altitude + # temperature + temperature = args.temperature + # Whether to use altitude + use_altitude = args.use_altitude + # Reynold's Length + reynolds_length = 1.0 + # volume constraint ratio + area_ratio_min = args.area_ratio_min + area_initial = args.area_initial + area_input_design = args.area_input_design + + # Optimization parameters + opt = args.opt + opt_options = args.opt_options + # ====================================================================== + # Create multipoint communication object + # ====================================================================== + mp = multiPointSparse(MPI.COMM_WORLD) + mp.addProcessorSet("cruise", nMembers=1, memberSizes=MPI.COMM_WORLD.size) + comm, *_ = mp.createCommunicators() + if not os.path.exists(args.output_dir) and comm.rank == 0: + os.mkdir(args.output_dir) + + # ====================================================================== + # ADflow Set-up + # ====================================================================== + aero_options = { + # Common Parameters + "gridFile": args.mesh_fname, + "outputDirectory": args.output_dir, + "writeVolumeSolution": False, + "writeTecplotSurfaceSolution": True, + "monitorvariables": ["cl", "cd", "yplus"], + # Physics Parameters + "equationType": "RANS", + "smoother": "DADI", + "nCycles": 5000, + "rkreset": True, + "nrkreset": 10, + # NK Options + "useNKSolver": True, + "nkswitchtol": 1e-8, + # ANK Options + "useanksolver": True, + "ankswitchtol": 1e-1, + "liftIndex": 2, + "infchangecorrection": True, + "nsubiterturb": 10, + # Convergence Parameters + "L2Convergence": 1e-8, + "L2ConvergenceCoarse": 1e-4, + # Adjoint Parameters + "adjointSolver": "GMRES", + "adjointL2Convergence": 1e-8, + "ADPC": True, + "adjointMaxIter": 1000, + "adjointSubspaceSize": 200, + } + + # Create solver + cfd_solver = ADFLOW(options=aero_options, comm=comm) + # ====================================================================== + # Set up flow conditions with AeroProblem + # ====================================================================== + + if use_altitude: + ap = AeroProblem( + name="fc", + alpha=alpha, + mach=mach, + altitude=altitude, + areaRef=1.0, + chordRef=1.0, + evalFuncs=["cl", "cd"], + ) + else: + ap = AeroProblem( + name="fc", + alpha=alpha, + mach=mach, + T=temperature, + reynolds=reynolds, + reynoldsLength=reynolds_length, + areaRef=1.0, + chordRef=1.0, + evalFuncs=["cl", "cd"], + ) + + # Add angle of attack variable + ap.addDV("alpha", value=alpha, lower=0.0, upper=10.0, scale=1.0) + # ====================================================================== + # Geometric Design Variable Set-up + # ====================================================================== + # Create DVGeometry object + dv_geo = DVGeometry(args.ffd_fname) + dv_geo.addLocalDV("shape", lower=-0.025, upper=0.025, axis="y", scale=1.0) + + span = 1.0 + pos = np.array([0.5]) * span + cfd_solver.addSlices("z", pos, sliceType="absolute") + + # Add DVGeo object to CFD solver + cfd_solver.setDVGeo(dv_geo) + # ====================================================================== + # DVConstraint Setup + # ====================================================================== + + dv_con = DVConstraints() + dv_con.setDVGeo(dv_geo) + + # Only ADflow has the getTriangulatedSurface Function + dv_con.setSurface(cfd_solver.getTriangulatedMeshSurface()) + + # Le/Te constraints + l_index = dv_geo.getLocalIndex(0) + ind_set_a = [] + ind_set_b = [] + for k in range(1): + ind_set_a.append(l_index[0, 0, k]) # all DV for upper and lower should be same but different sign + ind_set_b.append(l_index[0, 1, k]) + for k in range(1): + ind_set_a.append(l_index[-1, 0, k]) + ind_set_b.append(l_index[-1, 1, k]) + dv_con.addLeTeConstraints(0, indSetA=ind_set_a, indSetB=ind_set_b) + + # DV should be same along spanwise + l_index = dv_geo.getLocalIndex(0) + ind_set_a = [] + ind_set_b = [] + for i in range(l_index.shape[0]): + ind_set_a.append(l_index[i, 0, 0]) + ind_set_b.append(l_index[i, 0, 1]) + for i in range(l_index.shape[0]): + ind_set_a.append(l_index[i, 1, 0]) + ind_set_b.append(l_index[i, 1, 1]) + dv_con.addLinearConstraintsShape(ind_set_a, ind_set_b, factorA=1.0, factorB=-1.0, lower=0, upper=0) + + le = 0.010001 + le_list = [[le, 0, le], [le, 0, 1.0 - le]] + te_list = [[1.0 - le, 0, le], [1.0 - le, 0, 1.0 - le]] + + dv_con.addVolumeConstraint( + le_list, + te_list, + 2, + 100, + lower=area_ratio_min * area_initial / area_input_design, + upper=1.2 * area_initial / area_input_design, + scaled=True, + ) + dv_con.addThicknessConstraints2D(le_list, te_list, 2, 100, lower=0.15, upper=3.0) + # Final constraint to keep TE thickness at original or greater + dv_con.addThicknessConstraints1D(ptList=te_list, nCon=2, axis=[0, 1, 0], lower=1.0, scaled=True) + + if comm.rank == 0: + file_name = os.path.join(args.output_dir, "constraints.dat") + dv_con.writeTecplot(file_name) + # ====================================================================== + # Mesh Warping Set-up + # ====================================================================== + mesh_options = {"gridFile": args.mesh_fname} + + mesh = USMesh(options=mesh_options, comm=comm) + cfd_solver.setMesh(mesh) + + # ====================================================================== + # Optimization Problem Set-up + # ====================================================================== + # Create optimization problem + opt_prob = Optimization("opt", mp.obj, comm=MPI.COMM_WORLD) + + # Add objective + opt_prob.addObj("obj", scale=1e4) + + # Add variables from the AeroProblem + ap.addVariablesPyOpt(opt_prob) + + # Add DVGeo variables + dv_geo.addVariablesPyOpt(opt_prob) + + # Add constraints + dv_con.addConstraintsPyOpt(opt_prob) + opt_prob.addCon("cl_con_" + ap.name, lower=0.0, upper=0.0, scale=1.0) + + # The MP object needs the 'obj' and 'sens' function for each proc set, + # the optimization problem and what the objcon function is: + def cruise_funcs(x): + if MPI.COMM_WORLD.rank == 0: + print(x) + # Set design vars + dv_geo.setDesignVars(x) + ap.setDesignVars(x) + # Run CFD + cfd_solver(ap) + # Evaluate functions + funcs = {} + dv_con.evalFunctions(funcs) + cfd_solver.evalFunctions(ap, funcs) + cfd_solver.checkSolutionFailure(ap, funcs) + if MPI.COMM_WORLD.rank == 0: + print(funcs) + return funcs + + def cruise_funcs_sens(_x, _funcs): + funcs_sens = {} + dv_con.evalFunctionsSens(funcs_sens) + cfd_solver.evalFunctionsSens(ap, funcs_sens) + cfd_solver.checkAdjointFailure(ap, funcs_sens) + if MPI.COMM_WORLD.rank == 0: + print(funcs_sens) + return funcs_sens + + def obj_con(funcs, print_ok): + # Assemble the objective and any additional constraints: + funcs["obj"] = funcs[ap["cd"]] + funcs["cl_con_" + ap.name] = funcs[ap["cl"]] - mycl + if print_ok: + print("funcs in obj:", funcs) + return funcs + + mp.setProcSetObjFunc("cruise", cruise_funcs) + mp.setProcSetSensFunc("cruise", cruise_funcs_sens) + mp.setObjCon(obj_con) + mp.setOptProb(opt_prob) + opt_prob.printSparsity() + # Set up optimizer + if opt == Algorithm.SLSQP: + opt_options = {"IFILE": os.path.join(args.output_dir, "SLSQP.out")} + elif opt == Algorithm.SNOPT: + opt_options = { + "Major feasibility tolerance": 1e-4, + "Major optimality tolerance": 1e-4, + "Hessian full memory": None, + "Function precision": 1e-8, + "Print file": os.path.join(args.output_dir, "SNOPT_print.out"), + "Summary file": os.path.join(args.output_dir, "SNOPT_summary.out"), + } + opt_options.update(opt_options) + opt = OPT(opt.name, options=opt_options) + + # Run Optimization + sol = opt(opt_prob, mp.sens, sensMode="pgc", sensStep=1e-6, storeHistory=os.path.join(args.output_dir, "opt.hst")) + if MPI.COMM_WORLD.rank == 0: + print(sol) + + MPI.COMM_WORLD.Barrier() + + +if __name__ == "__main__": + main() diff --git a/engibench/problems/wings3D/templates/cli_interface.py b/engibench/problems/wings3D/templates/cli_interface.py new file mode 100644 index 00000000..b4768ece --- /dev/null +++ b/engibench/problems/wings3D/templates/cli_interface.py @@ -0,0 +1,129 @@ +"""Dataclasses sent from the main script to scripts inside the airfoil container.""" + +from dataclasses import asdict +from dataclasses import dataclass +from enum import auto +from enum import Enum +from typing import Any + + +class Task(Enum): + """The task to perform by analysis.""" + + ANALYSIS = auto() + POLAR = auto() + + +@dataclass +class AnalysisParameters: + """Parameters for airfoil_analyse.py.""" + + mesh_fname: str + """Path to the mesh file""" + output_dir: str + """Path to the output directory""" + task: Task + """The task to perform: "analysis" or "polar""" + mach: float + """mach number""" + reynolds: float + """Reynolds number""" + altitude: float + """altitude""" + temperature: float + """temperature""" + use_altitude: bool + """Whether to use altitude""" + alpha: float + """Angle of attack""" + + def to_dict(self) -> dict[str, Any]: + """Serialize to python primitives.""" + return {**asdict(self), "task": self.task.name} + + @classmethod + def from_dict(cls, data: dict[str, Any]) -> "AnalysisParameters": + """Deserialize from python primitives.""" + return cls(task=Task[data.pop("task")], **data) + + +class Algorithm(Enum): + """Algorithm to be used by optimize.""" + + SLSQP = auto() + SNOPT = auto() + + +@dataclass +class OptimizeParameters: + """Parameters for airfoil_opt.py.""" + + cl_target: float + """The lift coefficient constraint""" + alpha: float + """The angle of attack""" + mach: float + """The Mach number""" + reynolds: float + """Reynolds number""" + temperature: float + """Temperature""" + altitude: int + """The cruising altitude""" + area_ratio_min: float + area_initial: float + area_input_design: float + ffd_fname: str + """Path to the FFD file""" + mesh_fname: str + """Path to the mesh file""" + output_dir: str + """Path to the output directory""" + opt: Algorithm + """The optimization algorithm: SLSQP or SNOPT""" + opt_options: dict[str, Any] + """The optimization options""" + use_altitude: bool + """Whether to use altitude""" + + def to_dict(self) -> dict[str, Any]: + """Serialize to python primitives.""" + return {**asdict(self), "opt": self.opt.name} + + @classmethod + def from_dict(cls, data: dict[str, Any]) -> "OptimizeParameters": + """Deserialize from python primitives.""" + return cls(opt=Algorithm[data.pop("opt")], **data) + + +@dataclass +class PreprocessParameters: + """Parameters for pre_process.py.""" + + design_fname: str + """Path to the design file""" + N_sample: int + """Number of points to sample on the airfoil surface. Defines part of the mesh resolution""" + n_tept_s: int + """Number of points on the trailing edge""" + x_cut: float + """Blunt edge dimensionless cut location""" + tmp_xyz_fname: str + """Path to the temporary xyz file""" + mesh_fname: str + """Path to the generated mesh file""" + ffd_fname: str + """Path to the generated FFD file""" + ffd_ymarginu: float + """Upper (y-axis) margin for the fitted FFD cage""" + ffd_ymarginl: float + """Lower (y-axis) margin for the fitted FFD cage""" + ffd_pts: int + """Number of FFD points""" + N_grid: int + """Number of grid levels to march from the airfoil surface. Defines part of the mesh resolution""" + s0: float + """Off-the-wall spacing for the purpose of modeling the boundary layer""" + input_blunted: bool + march_dist: float + """Distance to march the grid from the airfoil surface""" diff --git a/engibench/problems/wings3D/templates/pre_process.py b/engibench/problems/wings3D/templates/pre_process.py new file mode 100644 index 00000000..ae0a978c --- /dev/null +++ b/engibench/problems/wings3D/templates/pre_process.py @@ -0,0 +1,75 @@ +# mypy: ignore-errors +"""This file is largely based on the MACHAero tutorials. + +https://github.com/mdolab/MACH-Aero/blob/main/tutorial/ + +TODO: Add the automatic grid spacing calculation. +""" + +import json +import sys + +from cli_interface import PreprocessParameters +import prefoil +from pyhyp import pyHyp + +if __name__ == "__main__": + args = PreprocessParameters(**json.loads(sys.argv[1])) + + coords = prefoil.utils.readCoordFile(args.design_fname) + airfoil = prefoil.Airfoil(coords) + print("Running pre-process.py") + input_blunted = args.input_blunted + if not input_blunted: + airfoil.normalizeAirfoil() + airfoil.makeBluntTE(xCut=args.x_cut) + + N_sample = args.N_sample + n_tept_s = args.n_tept_s + + coords = airfoil.getSampledPts( + N_sample, + spacingFunc=prefoil.sampling.conical, + func_args={"coeff": 1.2}, + nTEPts=n_tept_s, + ) + + # Write a fitted FFD with 10 chordwise points + ffd_ymarginu = args.ffd_ymarginu + ffd_ymarginl = args.ffd_ymarginl + ffd_fname = args.ffd_fname + ffd_pts = args.ffd_pts + airfoil.generateFFD(ffd_pts, ffd_fname, ymarginu=ffd_ymarginu, ymarginl=ffd_ymarginl) + + # write out plot3d + airfoil.writeCoords(args.tmp_xyz_fname, file_format="plot3d") + + # GenOptions + options = { + # --------------------------- + # Input Parameters + # --------------------------- + "inputFile": args.tmp_xyz_fname + ".xyz", + "unattachedEdgesAreSymmetry": False, + "outerFaceBC": "farfield", + "autoConnect": True, + "BC": {1: {"jLow": "zSymm", "jHigh": "zSymm"}}, + "families": "wall", + # --------------------------- + # Grid Parameters + # --------------------------- + "N": args.N_grid, + "nConstantStart": 8, + "s0": args.s0, + "marchDist": args.march_dist, + # Smoothing parameters + "volSmoothIter": 150, + "volCoef": 0.25, + "volBlend": 0.001, + } + + hyp = pyHyp(options=options) + hyp.run() + hyp.writeCGNS(args.mesh_fname) + + print(f"Generated files FFD and mesh in {ffd_fname}, {args.mesh_fname}") diff --git a/engibench/problems/wings3D/utils.py b/engibench/problems/wings3D/utils.py new file mode 100644 index 00000000..a1c6d2c6 --- /dev/null +++ b/engibench/problems/wings3D/utils.py @@ -0,0 +1,429 @@ +"""Utility functions for the airfoil problem.""" + +import numpy as np +import numpy.typing as npt +import pandas as pd + + +def _extract_connectivities(df_slice: pd.DataFrame) -> tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]: + """Extract node connectivities from the dataframe slice. + + Args: + df_slice (pd.DataFrame): A slice of a dataframe. + + Returns: + tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]: NodeC1 and NodeC2 arrays + """ + node_c1 = np.array(df_slice["NodeC1"].dropna().values).astype(int) # A list of node indices + node_c2 = np.array(df_slice["NodeC2"].dropna().values).astype(int) # A list of node indices + return node_c1, node_c2 + + +def _identify_segments(connectivities: npt.NDArray[np.int32]) -> tuple[list[int], list[int], npt.NDArray[np.float32]]: + """Identify segment breaks and assign segment IDs. + + Args: + connectivities (npt.NDArray[np.int32]): Array of node connections + + Returns: + tuple[list[int], list[int], npt.NDArray[np.float32]]: Start indices, end indices, and segment IDs + """ + id_breaks_start = [0] + id_breaks_end = [] + prev_id = 0 + segment_ids = np.zeros(len(connectivities), dtype=np.float32) + seg_id = 0 + j = -1 # Ensure j is always defined + + for j in range(len(connectivities)): + if connectivities[j][0] - 1 != prev_id: + # This means that we have a new set of points + id_breaks_start.append(connectivities[j][0] - 1) + id_breaks_end.append(prev_id) + seg_id += 1 + segment_ids[j] = seg_id + prev_id = connectivities[j][1] - 1 + + id_breaks_end.append(j) + return id_breaks_start, id_breaks_end, segment_ids + + +def _order_segments( + coords_x: npt.NDArray[np.float32], + coords_y: npt.NDArray[np.float32], + id_breaks_start: list[int], + id_breaks_end: list[int], + unique_segment_ids: npt.NDArray[np.int32], +) -> npt.NDArray[np.int32]: + """Order segments based on their spatial relationships. + + Args: + coords_x (npt.NDArray[np.float32]): X coordinates + coords_y (npt.NDArray[np.float32]): Y coordinates + id_breaks_start (list[int]): Start indices of segments + id_breaks_end (list[int]): End indices of segments + unique_segment_ids (npt.NDArray[np.int32]): Unique segment identifiers + + Returns: + npt.NDArray[np.int32]: Ordered segment IDs + """ + seg_coords_start_x = coords_x[id_breaks_start] + seg_coords_start_y = coords_y[id_breaks_start] + seg_coords_end_x = coords_x[id_breaks_end] + seg_coords_end_y = coords_y[id_breaks_end] + + ordered_ids = [unique_segment_ids[0]] + seg_coords_end_x_idx = seg_coords_end_x[0] + seg_coords_end_y_idx = seg_coords_end_y[0] + seg_coords_start_x_idx = seg_coords_start_x[0] + seg_coords_start_y_idx = seg_coords_start_y[0] + + # Loop through and find the end or start of a segment that matches the start of the current segment + while len(ordered_ids) < len(unique_segment_ids): + # Calculate the distance between the end of the current segment and the start of all other segments + diff_end_idx_start_tot = np.sqrt( + np.square(seg_coords_end_x_idx - seg_coords_start_x) + np.square(seg_coords_end_y_idx - seg_coords_start_y) + ) + diff_start_idx_start_tot = np.sqrt( + np.square(seg_coords_start_x_idx - seg_coords_start_x) + np.square(seg_coords_start_y_idx - seg_coords_start_y) + ) + + # Get the minimum distance excluding the current ordered segments + diff_end_idx_start_tot[np.abs(ordered_ids)] = np.inf + diff_start_idx_start_tot[np.abs(ordered_ids)] = np.inf + diff_end_idx_start_tot_id = np.argmin(diff_end_idx_start_tot) + diff_end_idx_start_tot_min = diff_end_idx_start_tot[diff_end_idx_start_tot_id] + + # Get the minimum distance excluding the current segment + diff_end_idx_end_tot = np.sqrt( + np.square(seg_coords_end_x_idx - seg_coords_end_x) + np.square(seg_coords_end_y_idx - seg_coords_end_y) + ) + diff_end_idx_end_tot[np.abs(ordered_ids)] = np.inf + diff_end_idx_end_tot_id = np.argmin(diff_end_idx_end_tot) + diff_end_idx_end_tot_min = diff_end_idx_end_tot[diff_end_idx_end_tot_id] + + # If the end of the current segment matches the start of another segment, + # we have found the correct order + if diff_end_idx_start_tot_min < diff_end_idx_end_tot_min: + # Append the matching segment id to the ordered ids + ordered_ids.append(diff_end_idx_start_tot_id) + # Update the current segment end coordinates + seg_coords_end_x_idx = seg_coords_end_x[diff_end_idx_start_tot_id] + seg_coords_end_y_idx = seg_coords_end_y[diff_end_idx_start_tot_id] + else: + # If the end of the current segment matches the end of another segment, + # the segment we append must be in reverse order + # We make the sign of the segment id negative to indicate reverse order + ordered_ids.append(-diff_end_idx_end_tot_id) + # Update the current segment end coordinates; + # Because of reversal, we use the start of the segment we are appending as the new end coordinates + seg_coords_end_x_idx = seg_coords_start_x[diff_end_idx_end_tot_id] + seg_coords_end_y_idx = seg_coords_start_y[diff_end_idx_end_tot_id] + + return np.array(ordered_ids) + + +def _reorder_coordinates( # noqa: PLR0913 + coords_x: npt.NDArray[np.float32], + coords_y: npt.NDArray[np.float32], + indices: npt.NDArray[np.int32], + connectivities: npt.NDArray[np.int32], + segment_ids: npt.NDArray[np.float32], + new_seg_order: npt.NDArray[np.int32], +) -> tuple[npt.NDArray[np.float32], npt.NDArray[np.float32], npt.NDArray[np.int32]]: + """Reorder coordinates based on segment order. + + Args: + coords_x (npt.NDArray[np.float32]): X coordinates + coords_y (npt.NDArray[np.float32]): Y coordinates + indices (npt.NDArray[np.int32]): Original indices + connectivities (npt.NDArray[np.int32]): Node connections + segment_ids (npt.NDArray[np.float32]): Segment identifiers + new_seg_order (npt.NDArray[np.int32]): New segment order + + Returns: + tuple[npt.NDArray[np.float32], npt.NDArray[np.float32], npt.NDArray[np.int32]]: Reordered coordinates and indices + """ + coords_x_reordered = np.array([]) + coords_y_reordered = np.array([]) + indices_reordered = np.array([]) + + for j in range(len(new_seg_order)): + if new_seg_order[j] < 0: + segment = np.nonzero(segment_ids == -new_seg_order[j])[0] + coords_x_segment = coords_x[connectivities[segment] - 1][:, 0][::-1] + coords_y_segment = coords_y[connectivities[segment] - 1][:, 0][::-1] + indices_segment = indices[connectivities[segment] - 1][:, 0][::-1] + else: + segment = np.nonzero(segment_ids == new_seg_order[j])[0] + coords_x_segment = coords_x[connectivities[segment] - 1][:, 0] + coords_y_segment = coords_y[connectivities[segment] - 1][:, 0] + indices_segment = indices[connectivities[segment] - 1][:, 0] + + coords_x_reordered = np.concatenate((coords_x_reordered, coords_x_segment)) + coords_y_reordered = np.concatenate((coords_y_reordered, coords_y_segment)) + indices_reordered = np.concatenate((indices_reordered, indices_segment)) + + return coords_x_reordered, coords_y_reordered, indices_reordered + + +def _align_coordinates( + coords_x_reordered: npt.NDArray[np.float32], + coords_y_reordered: npt.NDArray[np.float32], + indices_reordered: npt.NDArray[np.int32], + err: float = 1e-4, + err_x: float = 0.90, +) -> tuple[npt.NDArray[np.float32], npt.NDArray[np.float32], npt.NDArray[np.int32]]: + """Align coordinates based on maximum x values and mean y values. + + Args: + coords_x_reordered (npt.NDArray[np.float32]): Reordered x coordinates + coords_y_reordered (npt.NDArray[np.float32]): Reordered y coordinates + indices_reordered (npt.NDArray[np.int32]): Reordered indices + err (float): Error tolerance + err_x (float): X coordinate error factor + + Returns: + tuple[npt.NDArray[np.float32], npt.NDArray[np.float32], npt.NDArray[np.int32]]: Aligned coordinates and indices + """ + max_x = np.amax(coords_x_reordered) * err_x + max_x_ids = np.argwhere(np.abs(coords_x_reordered - np.amax(coords_x_reordered)) < err).reshape(-1, 1) + # Maintain (N,1) shape by using boolean indexing and reshaping + mask = coords_x_reordered[max_x_ids.ravel()] >= max_x + max_x_ids = max_x_ids[mask].reshape(-1, 1) + + # Get the y values at the maximum x values + max_x_y_values = coords_y_reordered[max_x_ids] + max_y_value = np.max(max_x_y_values) + min_y_value = np.min(max_x_y_values) + mean_y_value = (max_y_value + min_y_value) / 2 + + # Get the id of the value closest to the mean y value at the x value of the maximum y value + mean_y_value_sub_id = np.argmin(np.abs(max_x_y_values - mean_y_value)) + mean_y_value_id = max_x_ids[mean_y_value_sub_id].item() + + # Now reorder the coordinates such that the mean y value is first + coords_x_reordered = np.concatenate((coords_x_reordered[mean_y_value_id:], coords_x_reordered[:mean_y_value_id])) + coords_y_reordered = np.concatenate((coords_y_reordered[mean_y_value_id:], coords_y_reordered[:mean_y_value_id])) + indices_reordered = np.concatenate((indices_reordered[mean_y_value_id:], indices_reordered[:mean_y_value_id])) + + return coords_x_reordered, coords_y_reordered, indices_reordered + + +def _clean_coordinates( + coords_x_reordered: npt.NDArray[np.float32], + coords_y_reordered: npt.NDArray[np.float32], + indices_reordered: npt.NDArray[np.int32], + err_remove: float = 1e-6, +) -> tuple[npt.NDArray[np.float32], npt.NDArray[np.float32], npt.NDArray[np.int32]]: + """Remove duplicate coordinates and close the loop. + + Args: + coords_x_reordered (npt.NDArray[np.float32]): Reordered x coordinates + coords_y_reordered (npt.NDArray[np.float32]): Reordered y coordinates + indices_reordered (npt.NDArray[np.int32]): Reordered indices + err_remove (float): Error tolerance for removal + + Returns: + tuple[npt.NDArray[np.float32], npt.NDArray[np.float32], npt.NDArray[np.int32]]: Cleaned coordinates and indices + """ + removal_ids = np.where(np.abs(np.diff(coords_x_reordered) + np.diff(coords_y_reordered)) < err_remove)[0] + indices_reordered = np.delete(indices_reordered, removal_ids) + coords_x_reordered = np.delete(coords_x_reordered, removal_ids) + coords_y_reordered = np.delete(coords_y_reordered, removal_ids) + + coords_x_reordered = np.concatenate((coords_x_reordered, np.expand_dims(coords_x_reordered[0], axis=0))) + coords_y_reordered = np.concatenate((coords_y_reordered, np.expand_dims(coords_y_reordered[0], axis=0))) + indices_reordered = np.concatenate((indices_reordered, np.expand_dims(indices_reordered[0], axis=0))) + + return coords_x_reordered, coords_y_reordered, indices_reordered + + +def reorder_coords(df_slice: pd.DataFrame) -> npt.NDArray[np.float32]: + """Reorder the coordinates of a slice of a dataframe. + + Args: + df_slice (pd.DataFrame): A slice of a dataframe. + + Returns: + npt.NDArray[np.float32]: The reordered coordinates. + """ + # Extract connectivities + node_c1, node_c2 = _extract_connectivities(df_slice) + connectivities = np.concatenate((node_c1.reshape(-1, 1), node_c2.reshape(-1, 1)), axis=1) + + # Get coordinates + coords_x = df_slice["CoordinateX"].to_numpy() + coords_y = df_slice["CoordinateY"].to_numpy() + indices = np.arange(len(df_slice), dtype=np.int32) + + # Identify segments + id_breaks_start, id_breaks_end, segment_ids = _identify_segments(connectivities) + unique_segment_ids = np.arange(len(id_breaks_start), dtype=np.int32) + + # Order segments + new_seg_order = _order_segments(coords_x, coords_y, id_breaks_start, id_breaks_end, unique_segment_ids) + + # Reorder coordinates + coords_x_reordered, coords_y_reordered, indices_reordered = _reorder_coordinates( + coords_x, coords_y, indices, connectivities, segment_ids, new_seg_order + ) + + # Align coordinates + coords_x_reordered, coords_y_reordered, indices_reordered = _align_coordinates( + coords_x_reordered, coords_y_reordered, indices_reordered + ) + + # Clean coordinates + coords_x_reordered, coords_y_reordered, indices_reordered = _clean_coordinates( + coords_x_reordered, coords_y_reordered, indices_reordered + ) + + return np.array([coords_x_reordered, coords_y_reordered]) + + +def scale_coords( + coords: npt.NDArray[np.float64], + blunted: bool = False, # noqa: FBT001, FBT002 + xcut: float = 0.99, + min_trailing_edge_indices: float = 6, +) -> tuple[npt.NDArray[np.float64], bool]: + """Scales the coordinates to fit in the design space. + + Args: + coords (np.ndarray): The coordinates to scale. + blunted (bool): If True, the coordinates are assumed to be blunted. + xcut (float): The x coordinate of the cut, if the coordinates are blunted. + min_trailing_edge_indices (int): The minimum number of trailing edge indices to remove. + + Returns: + np.ndarray: The scaled coordinates. + """ + # Test if the coordinates are blunted or not + if not (blunted) and is_blunted(coords): + blunted = True + print( + "The coordinates may be blunted. However, blunted was not set to True. Will set blunted to True and continue, but please check the coordinates." + ) + + if not (blunted): + xcut = 1.0 + + # Scale x coordinates to be xcut in length + airfoil_length = np.abs(np.max(coords[0, :]) - np.min(coords[0, :])) + + # Center the coordinates around the leading edge and scale them + coords[0, :] = xcut * (coords[0, :] - np.min(coords[0, :])) / airfoil_length + airfoil_length = np.abs(np.max(coords[0, :]) - np.min(coords[0, :])) + + # Shift the coordinates to be centered at 0 at the leading edge + leading_id = np.argmin(coords[0, :]) + y_dist = coords[1, leading_id] + coords[1, :] += -y_dist + # Ensure the first and last points are the same + coords[0, 0] = xcut + coords[0, -1] = xcut + coords[1, -1] = coords[1, 0] + # Set the leading edge location + + if blunted: + coords_x = coords[0, :] + # Get all of the trailing edge indices, i.e where consecutive x coordinates are the same + err = 1e-5 + x_gt = np.max(coords_x) * 0.99 + trailing_edge_indices_l = np.where(np.abs(coords_x - np.roll(coords_x, -1)) < err)[0] + trailing_edge_indices_r = np.where(np.abs(coords_x - np.roll(coords_x, 1)) < err)[0] + # Include any indices that are in either list + trailing_edge_indices = np.unique(np.concatenate((trailing_edge_indices_l, trailing_edge_indices_r))) + trailing_edge_indices = trailing_edge_indices[coords_x[trailing_edge_indices] >= x_gt] + + err = 1e-4 + err_stop = 1e-3 + while len(trailing_edge_indices) < min_trailing_edge_indices: + trailing_edge_indices_l = np.where(np.abs(coords_x - np.roll(coords_x, -1)) < err)[0] + trailing_edge_indices_r = np.where(np.abs(coords_x - np.roll(coords_x, 1)) < err)[0] + # Include any indices that are in either list + trailing_edge_indices = np.unique(np.concatenate((trailing_edge_indices_l, trailing_edge_indices_r))) + trailing_edge_indices = trailing_edge_indices[coords_x[trailing_edge_indices] >= x_gt] + err *= 1.5 + if err > err_stop: + break + + # Remove the trailing edge indices from the coordinates + coords = np.delete(coords, trailing_edge_indices[1:-1], axis=1) + + return coords, blunted + + +def calc_off_wall_distance( # noqa: PLR0913 + mach: float, + reynolds: float, + freestreamTemp: float = 300.0, # noqa: N803 + reynoldsLength: float = 1.0, # noqa: N803 + yplus: float = 1, + R: float = 287.0, # noqa: N803 + gamma: float = 1.4, +) -> float: + """Estimation of the off-wall distance for a given design. + + The off-wall distance is calculated using the Reynolds number and the freestream temperature. + """ + # --------------------------- + a = np.sqrt(gamma * R * freestreamTemp) + u = mach * a + # --------------------------- + # Viscosity from Sutherland's law + ## Sutherland's law parameters + mu0 = 1.716e-5 + T0 = 273.15 # noqa: N806 + S = 110.4 # noqa: N806 + mu = mu0 * ((freestreamTemp / T0) ** (3 / 2)) * (T0 + S) / (freestreamTemp + S) + # --------------------------- + # Density + rho = reynolds * mu / (reynoldsLength * u) + ## Skin friction coefficient + Cf = (2 * np.log10(reynolds) - 0.65) ** (-2.3) # noqa: N806 + # Wall shear stress + tau = Cf * 0.5 * rho * (u**2) + # Friction velocity + uTau = np.sqrt(tau / rho) # noqa: N806 + # Off wall distance + return yplus * mu / (rho * uTau) + + +def is_blunted(coords: npt.NDArray[np.float64], delta_x_tol: float = 1e-5) -> bool: + """Checks if the coordinates are blunted or not. + + Args: + coords (np.ndarray): The coordinates to check. + delta_x_tol (float): The tolerance for the x coordinate difference. + + Returns: + bool: True if the coordinates are blunted, False otherwise. + """ + # Check if the coordinates going away from the tip have a small delta y + coords_x = coords[0, :] + # Get all of the trailing edge indices, i.e where consecutive x coordinates are the same + x_gt = np.max(coords_x) * 0.99 + trailing_edge_indices_l = np.where(np.abs(coords_x - np.roll(coords_x, -1)) < delta_x_tol)[0] + trailing_edge_indices_r = np.where(np.abs(coords_x - np.roll(coords_x, 1)) < delta_x_tol)[0] + # Include any indices that are in either list + trailing_edge_indices = np.unique(np.concatenate((trailing_edge_indices_l, trailing_edge_indices_r))) + trailing_edge_indices = trailing_edge_indices[coords_x[trailing_edge_indices] >= x_gt] + + # check if we have no trailing edge indices + return not len(trailing_edge_indices) <= 1 + + +def calc_area(coords: npt.NDArray[np.float32]) -> float: + """Calculates the area of the airfoil. + + Args: + coords (np.ndarray): The coordinates of the airfoil. + + Returns: + float: The area of the airfoil. + """ + return 0.5 * np.absolute( + np.sum(coords[0, :] * np.roll(coords[1, :], -1)) - np.sum(coords[1, :] * np.roll(coords[0, :], -1)) + ) diff --git a/engibench/problems/wings3D/v0.py b/engibench/problems/wings3D/v0.py new file mode 100644 index 00000000..303270bd --- /dev/null +++ b/engibench/problems/wings3D/v0.py @@ -0,0 +1,534 @@ +"""Airfoil problem. + +Filename convention is that folder paths do not end with /. For example, /path/to/folder is correct, but /path/to/folder/ is not. +""" + +import dataclasses +from dataclasses import dataclass +from dataclasses import field +from importlib.util import find_spec +import json +import os +import shutil +import sys +from typing import Annotated, Any + +from gymnasium import spaces +from matplotlib.figure import Figure +import matplotlib.pyplot as plt +import numpy as np +import numpy.typing as npt +import pandas as pd + +from engibench.constraint import bounded +from engibench.constraint import constraint +from engibench.constraint import IMPL +from engibench.constraint import THEORY +from engibench.core import ObjectiveDirection +from engibench.core import OptiStep +from engibench.core import Problem +from engibench.problems.airfoil.pyopt_history import History +from engibench.problems.airfoil.templates import cli_interface +from engibench.problems.airfoil.utils import calc_area +from engibench.problems.airfoil.utils import calc_off_wall_distance +from engibench.problems.airfoil.utils import reorder_coords +from engibench.problems.airfoil.utils import scale_coords +from engibench.utils import container +from engibench.utils.files import clone_dir + +# Allow loading pyoptsparse histories even if pyoptsparse is not installed: +if find_spec("pyoptsparse") is None: + from engibench.problems.airfoil import fake_pyoptsparse + + sys.modules["pyoptsparse"] = fake_pyoptsparse + +DesignType = dict[str, Any] + + +def self_intersect(curve: npt.NDArray[np.float64]) -> tuple[int, npt.NDArray[np.float64], npt.NDArray[np.float64]] | None: + """Determines if two segments a and b intersect.""" + # intersection: find t such that (p + t dp - q) x dq = 0 with 0 <= t <= 1 + # and (q + s dq - p) x dp = 0, 0 <= s <= 1 + # dp x dq = 0 => parallel => no intersection + # + # t = (q-p) x dq / dp x dq + # s = (q-p) x dp / dp x dq + # + # Also use the fact that 2 consecutive segments always intersect (at their common point) + # => never check consecutive segments + segments = curve[1:] - curve[:-1] + n = segments.shape[0] + for i in range(n - 1): + p, dp = curve[i], segments[i] + end = n - 1 if i == 0 else n + q, dq = curve[i + 2 : end], segments[i + 2 : end] + x = np.cross(dp, dq) + parallel = x == 0.0 + t = np.cross(q[~parallel] - p, dq[~parallel]) / x[~parallel] + s = np.cross(q[~parallel] - p, dp) / x[~parallel] + if np.any((t >= 0.0) & (t <= 1.0) & (s >= 0.0) & (s <= 1.0)): + return i, p, curve[i + 1] + return None + + +@constraint(categories=IMPL) +def does_not_self_intersect(design: DesignType) -> None: + """Check if a curve has no self intersections.""" + intersection = self_intersect(design["coords"]) + assert intersection is None, ( + f"design: Curve does self intersect at segment {intersection[0]}: {intersection[1]} -- {intersection[2]}" + ) + + +class Airfoil(Problem[DesignType]): + r"""Airfoil 2D shape optimization problem. + + This problem simulates the performance of an airfoil in a 2D environment. An airfoil is represented by a set of 192 points that define its shape. The performance is evaluated by the [MACH-Aero](https://mdolab-mach-aero.readthedocs-hosted.com/en/latest/) simulator that computes the lift and drag coefficients of the airfoil. + """ + + version = 0 + objectives: tuple[tuple[str, ObjectiveDirection], ...] = ( + ("cd", ObjectiveDirection.MINIMIZE), + ("cl", ObjectiveDirection.MAXIMIZE), + ) + + design_space = spaces.Dict( + { + "coords": spaces.Box(low=0.0, high=1.0, shape=(2, 192), dtype=np.float32), + "angle_of_attack": spaces.Box(low=0.0, high=10.0, shape=(1,), dtype=np.float32), + } + ) + design_constraints = (does_not_self_intersect,) + dataset_id = "IDEALLab/airfoil_v0" + container_id = "mdolab/public:u22-gcc-ompi-stable" + __local_study_dir: str + + @dataclass + class Conditions: + """Conditions.""" + + mach: Annotated[ + float, bounded(lower=0.0).category(IMPL), bounded(lower=0.1, upper=1.0).warning().category(IMPL) + ] = 0.8 + """Mach number""" + reynolds: Annotated[ + float, bounded(lower=0.0).category(IMPL), bounded(lower=1e5, upper=1e9).warning().category(IMPL) + ] = 1e6 + """Reynolds number""" + area_initial: float = float("NAN") + """actual initial airfoil area""" + area_ratio_min: Annotated[float, bounded(lower=0.0, upper=1.2).category(THEORY)] = 0.7 + """Minimum ratio the initial area is allowed to decrease to i.e minimum_area = area_initial*area_target""" + cl_target: float = 0.5 + """Target lift coefficient to satisfy equality constraint""" + + conditions = Conditions() + + @dataclass + class Config(Conditions): + """Structured representation of configuration parameters for a numerical computation.""" + + alpha: Annotated[float, bounded(lower=0.0, upper=10.0).category(THEORY)] = 0.0 + altitude: float = 10000.0 + temperature: float = 300.0 + use_altitude: bool = False + output_dir: str | None = None + mesh_fname: str | None = None + task: str = "analysis" + opt: str = "SLSQP" + opt_options: dict = field(default_factory=dict) + ffd_fname: str | None = None + area_input_design: float | None = None + + @constraint(categories=THEORY) + @staticmethod + def area_ratio_bound(area_ratio_min: float, area_initial: float, area_input_design: float | None) -> None: + """Constraint for area_ratio_min <= area_ratio <= 1.2.""" + area_ratio_max = 1.2 + if area_input_design is None: + return + assert not np.isnan(area_initial) + area_ratio = area_input_design / area_initial + assert area_ratio_min <= area_ratio <= area_ratio_max, ( + f"Config.area_ratio: {area_ratio} ∉ [area_ratio_min={area_ratio_min}, 1.2]" + ) + + def __init__(self, seed: int = 0, base_directory: str | None = None) -> None: + """Initializes the Airfoil problem. + + Args: + seed (int): The random seed for the problem. + base_directory (str, optional): The base directory for the problem. If None, the current directory is selected. + """ + # This is used for intermediate files + # Local file are prefixed with self.local_base_directory + if base_directory is not None: + self.__local_base_directory = base_directory + else: + self.__local_base_directory = os.getcwd() + self.__local_target_dir = self.__local_base_directory + "/engibench_studies/problems/airfoil" + self.__local_template_dir = ( + os.path.dirname(os.path.abspath(__file__)) + "/templates" + ) # These templates are shipped with the lib + self.__local_scripts_dir = os.path.dirname(os.path.abspath(__file__)) + "/scripts" + + # Docker target directory + # This is used for files that are mounted into the docker container + self.__docker_base_dir = "/home/mdolabuser/mount/engibench" + self.__docker_target_dir = self.__docker_base_dir + "/engibench_studies/problems/airfoil" + + super().__init__(seed=seed) + + def reset(self, seed: int | None = None, *, cleanup: bool = False) -> None: + """Resets the simulator and numpy random to a given seed. + + Args: + seed (int, optional): The seed to reset to. If None, a random seed is used. + cleanup (bool): Deletes the previous study directory if True. + """ + if cleanup: + shutil.rmtree(self.__local_study_dir) + + super().reset(seed) + self.current_study = f"study_{self.seed}-pid{os.getpid()}" + self.__local_study_dir = self.__local_target_dir + "/" + self.current_study + self.__docker_study_dir = self.__docker_target_dir + "/" + self.current_study + + def __design_to_simulator_input( + self, design: DesignType, mach: float, reynolds: float, temperature: float, filename: str = "design" + ) -> str: + """Converts a design to a simulator input. + + The simulator inputs are two files: a mesh file (.cgns) and a FFD file (.xyz). This function generates these files from the design. + The files are saved in the current directory with the name "$filename.cgns" and "$filename_ffd.xyz". + + Args: + design (dict): The design to convert. + mach: mach number + reynolds: reynolds number + temperature: temperature + filename (str): The filename to save the design to. + """ + # Creates the study directory + clone_dir(source_dir=self.__local_template_dir, target_dir=self.__local_study_dir) + + tmp = os.path.join(self.__docker_study_dir, "tmp") + + # Calculate the off-the-wall distance + estimate_s0 = True + + s0 = calc_off_wall_distance(mach=mach, reynolds=reynolds, freestreamTemp=temperature) if estimate_s0 else 1e-5 + # Scale the design to fit in the design space + x_cut = 0.99 + scaled_design, input_blunted = scale_coords( + design["coords"], + blunted=False, + xcut=x_cut, + ) + args = cli_interface.PreprocessParameters( + design_fname=f"{self.__docker_study_dir}/{filename}.dat", + tmp_xyz_fname=tmp, + mesh_fname=self.__docker_study_dir + "/" + filename + ".cgns", + ffd_fname=self.__docker_study_dir + "/" + filename + "_ffd", + N_sample=180, + n_tept_s=4, + x_cut=x_cut, + ffd_ymarginu=0.05, + ffd_ymarginl=0.05, + ffd_pts=10, + N_grid=100, + s0=s0, + input_blunted=input_blunted, + march_dist=100.0, + ) + + # Save the design to a temporary file. Format to 1e-6 rounding + np.savetxt(self.__local_study_dir + "/" + filename + ".dat", scaled_design.transpose()) + + # Launches a docker container with the pre_process.py script + # The script generates the mesh and FFD files + bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && python {self.__docker_study_dir}/pre_process.py '{json.dumps(dataclasses.asdict(args))}'" + assert self.container_id is not None, "Container ID is not set" + container.run( + command=["/bin/bash", "-c", bash_command], + image=self.container_id, + name="machaero", + mounts=[(self.__local_base_directory, self.__docker_base_dir)], + sync_uid=True, + ) + + return filename + + def simulator_output_to_design(self, simulator_output: str | None = None) -> npt.NDArray[np.float32]: + """Converts a simulator output to a design. + + Args: + simulator_output (str): The simulator output to convert. If None, the latest slice file is used. + + Returns: + np.ndarray: The corresponding design. + """ + if simulator_output is None: + # Take latest slice file + files = os.listdir(self.__local_study_dir + "/output") + files = [f for f in files if f.endswith("_slices.dat")] + file_numbers = [int(f.split("_")[1]) for f in files] + simulator_output = files[file_numbers.index(max(file_numbers))] + + slice_file = self.__local_study_dir + "/output/" + simulator_output + + # Define the variable names for columns + var_names = [ + "CoordinateX", + "CoordinateY", + "CoordinateZ", + "XoC", + "YoC", + "ZoC", + "VelocityX", + "VelocityY", + "VelocityZ", + "CoefPressure", + "Mach", + ] + + nelems = pd.read_csv( + slice_file, sep=r"\s+", names=["fill1", "Nodes", "fill2", "Elements", "ZONETYPE"], skiprows=3, nrows=1 + ) + nnodes = int(nelems["Nodes"].iloc[0]) + + # Read the main data and node connections + slice_df = pd.read_csv(slice_file, sep=r"\s+", names=var_names, skiprows=5, nrows=nnodes, engine="c") + nodes_arr = pd.read_csv(slice_file, sep=r"\s+", names=["NodeC1", "NodeC2"], skiprows=5 + nnodes, engine="c") + + # Concatenate node connections to the main data + slice_df = pd.concat([slice_df, nodes_arr], axis=1) + + return reorder_coords(slice_df) + + def simulate(self, design: DesignType, config: dict[str, Any] | None = None, mpicores: int = 4) -> npt.NDArray: + """Simulates the performance of an airfoil design. + + Args: + design (dict): The design to simulate. + config (dict): A dictionary with configuration (e.g., boundary conditions, filenames) for the simulation. + mpicores (int): The number of MPI cores to use in the simulation. + + Returns: + dict: The performance of the design - each entry of the dict corresponds to a named objective value. + """ + if isinstance(design["angle_of_attack"], np.ndarray): + design["angle_of_attack"] = design["angle_of_attack"][0] + + # pre-process the design and run the simulation + + # Prepares the airfoil_analysis.py script with the simulation configuration + conditions = self.Conditions() + config = config or {} + args = cli_interface.AnalysisParameters( + alpha=design["angle_of_attack"], + altitude=config.get("altitude", 10000), + temperature=config.get("temperature", 300), + reynolds=config.get("reynolds", conditions.reynolds), + mach=config.get("mach", conditions.mach), + use_altitude=config.get("use_altitude", False), + output_dir=config.get("output_dir", self.__docker_study_dir + "/output/"), + mesh_fname=config.get("mesh_fname", self.__docker_study_dir + "/design.cgns"), + task=cli_interface.Task[config["task"]] if "task" in config else cli_interface.Task.ANALYSIS, + ) + self.__design_to_simulator_input(design, mach=args.mach, reynolds=args.reynolds, temperature=args.temperature) + + # Launches a docker container with the airfoil_analysis.py script + # The script takes a mesh and ffd and performs an optimization + bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && mpirun -np {mpicores} python -m mpi4py {self.__docker_study_dir}/airfoil_analysis.py '{json.dumps(args.to_dict())}'" + assert self.container_id is not None, "Container ID is not set" + container.run( + command=["/bin/bash", "-c", bash_command], + image=self.container_id, + name="machaero", + mounts=[(self.__local_base_directory, self.__docker_base_dir)], + sync_uid=True, + ) + + outputs = np.load(self.__local_study_dir + "/output/outputs.npy") + lift = float(outputs[3]) + drag = float(outputs[4]) + return np.array([drag, lift]) + + def optimize( + self, starting_point: DesignType, config: dict[str, Any] | None = None, mpicores: int = 4 + ) -> tuple[DesignType, list[OptiStep]]: + """Optimizes the design of an airfoil. + + Args: + starting_point (dict): The starting point for the optimization. + config (dict): A dictionary with configuration (e.g., boundary conditions, filenames) for the optimization. + mpicores (int): The number of MPI cores to use in the optimization. + + Returns: + tuple[dict[str, Any], list[OptiStep]]: The optimized design and its performance. + """ + if isinstance(starting_point["angle_of_attack"], np.ndarray): + starting_point["angle_of_attack"] = starting_point["angle_of_attack"][0] + + # pre-process the design and run the simulation + filename = "candidate_design" + + # Prepares the optimize_airfoil.py script with the optimization configuration + fields = {f.name for f in dataclasses.fields(cli_interface.OptimizeParameters)} + config = {key: val for key, val in (config or {}).items() if key in fields} + if "area_initial" not in config: + raise ValueError("optimize(): config is missing the required parameter 'area_initial'") + if "opt" in config: + config["opt"] = cli_interface.Algorithm[config["opt"]] + args = cli_interface.OptimizeParameters( + **{ + **dataclasses.asdict(self.Conditions()), + "alpha": starting_point["angle_of_attack"], + "altitude": 10000, + "temperature": 300, # should specify either mach + altitude or mach + reynolds + reynoldsLength (default to 1) + temperature + "use_altitude": False, + "opt": cli_interface.Algorithm.SLSQP, + "opt_options": {}, + "output_dir": self.__docker_study_dir + "/output", + "ffd_fname": self.__docker_study_dir + "/" + filename + "_ffd.xyz", + "mesh_fname": self.__docker_study_dir + "/" + filename + ".cgns", + "area_input_design": calc_area(starting_point["coords"]), + **config, + }, + ) + self.__design_to_simulator_input( + starting_point, reynolds=args.reynolds, mach=args.reynolds, temperature=args.temperature, filename=filename + ) + + # Launches a docker container with the optimize_airfoil.py script + # The script takes a mesh and ffd and performs an optimization + bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && mpirun -np {mpicores} python -m mpi4py {self.__docker_study_dir}/airfoil_opt.py '{json.dumps(args.to_dict())}'" + assert self.container_id is not None, "Container ID is not set" + container.run( + command=["/bin/bash", "-c", bash_command], + image=self.container_id, + name="machaero", + mounts=[(self.__local_base_directory, self.__docker_base_dir)], + sync_uid=True, + ) + + # post process -- extract the shape and objective values + optisteps_history = [] + history = History(self.__local_study_dir + "/output/opt.hst") + call_counters = history.getCallCounters() + iters = list(map(int, call_counters)) if call_counters is not None else [] + + for i in range(len(iters)): + vals = history.read(int(iters[i])) + if vals is not None and "funcs" in vals and "obj" in vals["funcs"] and not vals["fail"]: + values = history.getValues(names=["obj"], callCounters=[i], allowSens=False, major=False, scale=True) + if values is not None and "obj" in values: + objective = values["obj"] + # flatten objective if it is a list + obj_np = np.array(objective) + if obj_np.ndim > 1: + obj_np = obj_np.flatten() + optisteps_history.append(OptiStep(obj_values=obj_np, step=vals["iter"])) + + history.close() + + opt_coords = self.simulator_output_to_design() + + return {"coords": opt_coords, "angle_of_attack": starting_point["angle_of_attack"]}, optisteps_history + + def render(self, design: DesignType, *, open_window: bool = False, save: bool = False) -> Figure: + """Renders the design in a human-readable format. + + Args: + design (dict): The design to render. + open_window (bool): If True, opens a window with the rendered design. + save (bool): If True, saves the rendered design to a file in the study directory. + + Returns: + Figure: The rendered design. + """ + fig, ax = plt.subplots() + coords = design["coords"] + alpha = design["angle_of_attack"] + ax.scatter(coords[0], coords[1], s=10, alpha=0.7) + ax.set_title(r"$\alpha$=" + str(np.round(alpha, 2)) + r"$^\circ$") + ax.axis("equal") + ax.axis("off") + ax.set_xlim((-0.005, 1.005)) + + if open_window: + plt.show() + if save: + plt.savefig(self.__local_study_dir + "/airfoil.png", dpi=300, bbox_inches="tight") + plt.close(fig) + return fig + + def render_optisteps(self, optisteps_history: list[OptiStep], *, open_window: bool = False, save: bool = False) -> Any: + """Renders the optimization step history. + + Args: + optisteps_history (list[OptiStep]): The optimization steps to render. + open_window (bool): If True, opens a window with the rendered design. + save (bool): If True, saves the rendered design to a file in the study directory. + + Returns: + Any: Rendered optimization step history. + """ + fig, ax = plt.subplots() + steps = np.array([step.step for step in optisteps_history]) + objectives = np.array([step.obj_values[0][0] for step in optisteps_history]) + ax.plot(steps, objectives, label="Drag Coefficient") + ax.set_title("Optimization Steps") + ax.set_xlabel("Iteration") + ax.set_ylabel("Drag counts") + if open_window: + plt.show() + if save: + plt.savefig(self.__local_study_dir + "/optisteps.png", dpi=300, bbox_inches="tight") + plt.close(fig) + return fig, ax + + def random_design(self, dataset_split: str = "train", design_key: str = "initial_design") -> tuple[dict[str, Any], int]: + """Samples a valid random initial design. + + Args: + dataset_split (str): The key to use for the dataset. Defaults to "train". + design_key (str): The key to use for the design in the dataset. + Defaults to "initial_design". + + Returns: + tuple[dict[str, Any], int]: The valid random design and the index of the design in the dataset. + """ + rnd = self.np_random.integers(low=0, high=len(self.dataset[dataset_split][design_key]), dtype=int) + initial_design = self.dataset[dataset_split][design_key][rnd] + return {"coords": np.array(initial_design["coords"]), "angle_of_attack": initial_design["angle_of_attack"]}, rnd + + +if __name__ == "__main__": + # Initialize the problem + + problem = Airfoil(seed=0) + + # Retrieve the dataset + dataset = problem.dataset + + # Get random initial design and optimized conditions from the dataset + the index + design, idx = problem.random_design() + + # Get the config conditions from the dataset + config = dataset["train"].select_columns(problem.conditions_keys)[idx] + + # Simulate the design + print("Simulation results: ", problem.simulate(design, config=config, mpicores=8)) + + # Cleanup the study directory; will delete the previous contents from simulate in this case + problem.reset(seed=1, cleanup=True) + + # Get design and conditions from the dataset, render design + opt_design, optisteps_history = problem.optimize(design, config=config, mpicores=8) + print("Optimized design: ", opt_design) + print("Optimization history: ", optisteps_history) + + # Render the final optimized design + problem.render(opt_design, open_window=False, save=True) From f98e9fd9a47f3eec877e023eb7a8aa47d9c47d2e Mon Sep 17 00:00:00 2001 From: Anna Delbeke Date: Wed, 4 Mar 2026 14:10:23 +0100 Subject: [PATCH 02/10] Make Wings3D HF loader robust --- .../problems/wings3D/dataset_hf_wings3d.py | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 engibench/problems/wings3D/dataset_hf_wings3d.py diff --git a/engibench/problems/wings3D/dataset_hf_wings3d.py b/engibench/problems/wings3D/dataset_hf_wings3d.py new file mode 100644 index 00000000..cdfeb4fd --- /dev/null +++ b/engibench/problems/wings3D/dataset_hf_wings3d.py @@ -0,0 +1,22 @@ +""" +HuggingFace dataset loader for the Wings3D problem. +""" + +from __future__ import annotations + +from datasets import DatasetDict, load_dataset + + +def load_wings3d_dataset(dataset_id: str) -> DatasetDict: + if not dataset_id: + raise ValueError("dataset_id must be a non-empty string") + + try: + return load_dataset(dataset_id) + except Exception as e: + raise RuntimeError( + f"Could not load Hugging Face dataset '{dataset_id}'.\n" + f"- If it hasn't been uploaded yet, this is expected.\n" + f"- If it's private, run: huggingface-cli login\n" + f"Original error: {type(e).__name__}: {e}" + ) from e \ No newline at end of file From 18978d3977d122ab53d65d958a048ef90f002bc6 Mon Sep 17 00:00:00 2001 From: Anna Delbeke Date: Wed, 4 Mar 2026 14:23:08 +0100 Subject: [PATCH 03/10] Add dataset-backed Wings3DOffline using slice extraction --- engibench/problems/wings3D/utils_wing3d.py | 0 engibench/problems/wings3D/v0_offline.py | 0 2 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 engibench/problems/wings3D/utils_wing3d.py create mode 100644 engibench/problems/wings3D/v0_offline.py diff --git a/engibench/problems/wings3D/utils_wing3d.py b/engibench/problems/wings3D/utils_wing3d.py new file mode 100644 index 00000000..e69de29b diff --git a/engibench/problems/wings3D/v0_offline.py b/engibench/problems/wings3D/v0_offline.py new file mode 100644 index 00000000..e69de29b From ca6fbaa9c8f41ea5efb64447657a1ad3e7233405 Mon Sep 17 00:00:00 2001 From: Anna Delbeke Date: Wed, 4 Mar 2026 14:24:54 +0100 Subject: [PATCH 04/10] Fix Wings3DOffline import path in problems package init --- engibench/problems/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/engibench/problems/__init__.py b/engibench/problems/__init__.py index 069a7ff8..9fa6938c 100644 --- a/engibench/problems/__init__.py +++ b/engibench/problems/__init__.py @@ -1 +1,2 @@ """Contains all the different problems modeled in the library.""" +from .wings3D.v0_offline import Wings3DOffline \ No newline at end of file From c5a323382e659846631a6c24758588522c03577a Mon Sep 17 00:00:00 2001 From: Anna Delbeke Date: Wed, 4 Mar 2026 14:28:44 +0100 Subject: [PATCH 05/10] Fix render return in Wings3DOffline --- engibench/problems/wings3D/v0_offline.py | 74 ++++++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/engibench/problems/wings3D/v0_offline.py b/engibench/problems/wings3D/v0_offline.py index e69de29b..0e1b0284 100644 --- a/engibench/problems/wings3D/v0_offline.py +++ b/engibench/problems/wings3D/v0_offline.py @@ -0,0 +1,74 @@ +""" +Wings3D (offline) dataset-backed problem. + +One sample corresponds to one slice (selected via slice_num) of a wing. +simulate() returns cd/cl stored in the dataset row. +""" + +from __future__ import annotations + +from typing import Any + +import numpy as np +import numpy.typing as npt +from gymnasium import spaces +from matplotlib.figure import Figure +import matplotlib.pyplot as plt + +from engibench.core import ObjectiveDirection, Problem +from engibench.problems.wings3D.dataset_hf_wings3d import load_wings3d_dataset +from engibench.problems.wings3D.utils_wing3d import get_slice_coords + +DesignType = dict[str, Any] + + +class Wings3DOffline(Problem[DesignType]): + version = 0 + objectives = (("cd", ObjectiveDirection.MINIMIZE), ("cl", ObjectiveDirection.MAXIMIZE)) + + design_space = spaces.Dict( + { + "coords": spaces.Box(low=-1e6, high=1e6, shape=(192, 2), dtype=np.float32), + "alpha": spaces.Box(low=-90.0, high=90.0, shape=(1,), dtype=np.float32), + } + ) + + dataset_id: str = "IDEALLab/wings3d_v0" # replace with real HF id when available + + def __init__(self, seed: int = 0) -> None: + super().__init__(seed=seed) + self.dataset = load_wings3d_dataset(self.dataset_id) + + def random_design(self, dataset_split: str = "train") -> tuple[DesignType, dict[str, Any], int]: + ds = self.dataset[dataset_split] + idx = int(self.np_random.integers(low=0, high=len(ds), dtype=int)) + row = ds[idx] + + slice_num = int(row["slice_num"]) + coords_slice = get_slice_coords(row["coords"], slice_num) + alpha = np.asarray([row.get("alpha", 0.0)], dtype=np.float32) + + design: DesignType = {"coords": coords_slice, "alpha": alpha} + config = {"split": dataset_split, "idx": idx} + return design, config, idx + + def simulate(self, design: DesignType, config: dict[str, Any] | None = None, **kwargs) -> npt.NDArray[np.float64]: + if not config or "split" not in config or "idx" not in config: + raise ValueError("simulate() requires config with keys {'split','idx'} (use random_design()).") + + row = self.dataset[config["split"]][int(config["idx"])] + cd = float(row.get("cd_val", np.nan)) + cl = float(row.get("cl_val", np.nan)) + return np.array([cd, cl], dtype=np.float64) + + def render(self, design: DesignType, *, open_window: bool = False, save: bool = False) -> Figure: + coords = np.asarray(design["coords"], dtype=np.float32) + fig, ax = plt.subplots() + ax.plot(coords[:, 0], coords[:, 1]) + ax.axis("equal") + ax.axis("off") + if open_window: + plt.show() + plt.close(fig) + return fig + \ No newline at end of file From 25880e789f9dd8f3ae51cb9c29180e318a87ea3b Mon Sep 17 00:00:00 2001 From: Anna Delbeke Date: Wed, 4 Mar 2026 14:30:10 +0100 Subject: [PATCH 06/10] Fix indentation in Wings3DOffline render --- engibench/problems/wings3D/v0_offline.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engibench/problems/wings3D/v0_offline.py b/engibench/problems/wings3D/v0_offline.py index 0e1b0284..e6964dec 100644 --- a/engibench/problems/wings3D/v0_offline.py +++ b/engibench/problems/wings3D/v0_offline.py @@ -61,7 +61,7 @@ def simulate(self, design: DesignType, config: dict[str, Any] | None = None, **k cl = float(row.get("cl_val", np.nan)) return np.array([cd, cl], dtype=np.float64) - def render(self, design: DesignType, *, open_window: bool = False, save: bool = False) -> Figure: + def render(self, design: DesignType, *, open_window: bool = False, save: bool = False) -> Figure: coords = np.asarray(design["coords"], dtype=np.float32) fig, ax = plt.subplots() ax.plot(coords[:, 0], coords[:, 1]) From 907c3a593cceba1aec91300d7bc61e30c96fc7be Mon Sep 17 00:00:00 2001 From: Anna Delbeke Date: Wed, 4 Mar 2026 14:32:37 +0100 Subject: [PATCH 07/10] Fix utils_wing3d get_slice_coords import --- engibench/problems/wings3D/utils_wing3d.py | 26 ++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/engibench/problems/wings3D/utils_wing3d.py b/engibench/problems/wings3D/utils_wing3d.py index e69de29b..ee4836ec 100644 --- a/engibench/problems/wings3D/utils_wing3d.py +++ b/engibench/problems/wings3D/utils_wing3d.py @@ -0,0 +1,26 @@ +""" +Utility functions for the Wings3D problem. +""" + +from __future__ import annotations + +import numpy as np +import numpy.typing as npt + + +def get_slice_coords(coords: npt.NDArray, slice_num: int) -> npt.NDArray[np.float32]: + """ + Extract one slice curve from coords. + + Expected shapes: + - (9, 192, 2): full wing sections, return coords[slice_num] -> (192, 2) + - (192, 2): already a single slice, return as-is + """ + arr = np.asarray(coords, dtype=np.float32) + + if arr.ndim == 3: + return arr[int(slice_num)] + if arr.ndim == 2: + return arr + + raise ValueError(f"Unexpected coords shape: {arr.shape}") \ No newline at end of file From 32242c59b3e13aae243f1d50afd8ef9a51925976 Mon Sep 17 00:00:00 2001 From: Anna Delbeke Date: Fri, 13 Mar 2026 13:23:42 +0100 Subject: [PATCH 08/10] feat: implement Wings3D 3D wing optimization problem (with dummy simulator/optimizer) - Refactor Airfoil class to Wings3D with 9-slice 3D design space (9x192x2) - Add dataset loader with HuggingFace fallback to local .pkl/.parquet files - Update utils to handle 3D coordinates and slice operations - Add new CFD simulation templates (wings3D_analysis.py, wings3D_opt.py) - Update module imports, problem registration, and remove obsolete files - Change dataset ID to "Cashen/optiwing3d_engibench" - Note: simulate() and optimize() use dummy implementations until 3D CFD is integrated --- engibench/problems/__init__.py | 2 +- engibench/problems/photonics2d/__init__.py | 2 +- .../problems/wings3D/dataset_hf_wings3d.py | 50 +++- ...irfoil_analysis.py => wings3D_analysis.py} | 0 .../{airfoil_opt.py => wings3D_opt.py} | 0 engibench/problems/wings3D/utils.py | 158 +++++++---- engibench/problems/wings3D/utils_wing3d.py | 26 -- engibench/problems/wings3D/v0.py | 245 +++++++++--------- engibench/problems/wings3D/v0_offline.py | 74 ------ engibench/utils/all_problems.py | 35 ++- 10 files changed, 292 insertions(+), 300 deletions(-) rename engibench/problems/wings3D/templates/{airfoil_analysis.py => wings3D_analysis.py} (100%) rename engibench/problems/wings3D/templates/{airfoil_opt.py => wings3D_opt.py} (100%) delete mode 100644 engibench/problems/wings3D/utils_wing3d.py delete mode 100644 engibench/problems/wings3D/v0_offline.py diff --git a/engibench/problems/__init__.py b/engibench/problems/__init__.py index 9fa6938c..c619c079 100644 --- a/engibench/problems/__init__.py +++ b/engibench/problems/__init__.py @@ -1,2 +1,2 @@ """Contains all the different problems modeled in the library.""" -from .wings3D.v0_offline import Wings3DOffline \ No newline at end of file +from .wings3D.v0 import Wings3D \ No newline at end of file diff --git a/engibench/problems/photonics2d/__init__.py b/engibench/problems/photonics2d/__init__.py index 6e0efe45..9544c975 100644 --- a/engibench/problems/photonics2d/__init__.py +++ b/engibench/problems/photonics2d/__init__.py @@ -1,5 +1,5 @@ """Photonics2D problem module.""" -from engibench.problems.photonics2d.v0 import Photonics2D +#from engibench.problems.photonics2d.v0 import Photonics2D __all__ = ["Photonics2D"] diff --git a/engibench/problems/wings3D/dataset_hf_wings3d.py b/engibench/problems/wings3D/dataset_hf_wings3d.py index cdfeb4fd..cbcb3d6e 100644 --- a/engibench/problems/wings3D/dataset_hf_wings3d.py +++ b/engibench/problems/wings3D/dataset_hf_wings3d.py @@ -1,22 +1,54 @@ """ -HuggingFace dataset loader for the Wings3D problem. +Dataset loader for the Wings3D problem. + +Tries Hugging Face first. Optionally falls back to a local file for development. """ from __future__ import annotations -from datasets import DatasetDict, load_dataset +from pathlib import Path +from typing import Optional + +import numpy as np +import pandas as pd +from datasets import Dataset, DatasetDict, load_dataset + + +def _pandas_to_datasetdict(df: pd.DataFrame, split: str = "train") -> DatasetDict: + # Convert numpy arrays to lists so HF Dataset can serialize them + df2 = df.copy() + if "coords" in df2.columns: + df2["coords"] = df2["coords"].apply(lambda x: np.asarray(x).tolist()) + return DatasetDict({split: Dataset.from_pandas(df2, preserve_index=False)}) -def load_wings3d_dataset(dataset_id: str) -> DatasetDict: +def load_wings3d_dataset(dataset_id: str, local_path: Optional[str] = None) -> DatasetDict: if not dataset_id: raise ValueError("dataset_id must be a non-empty string") try: return load_dataset(dataset_id) except Exception as e: - raise RuntimeError( - f"Could not load Hugging Face dataset '{dataset_id}'.\n" - f"- If it hasn't been uploaded yet, this is expected.\n" - f"- If it's private, run: huggingface-cli login\n" - f"Original error: {type(e).__name__}: {e}" - ) from e \ No newline at end of file + if local_path is None: + raise RuntimeError( + f"Could not load Hugging Face dataset '{dataset_id}'.\n" + f"- If it hasn't been uploaded yet, this is expected.\n" + f"- If it's private, run: huggingface-cli login\n" + f"Original error: {type(e).__name__}: {e}" + ) from e + + p = Path(local_path) + if not p.exists(): + raise RuntimeError( + f"HF dataset '{dataset_id}' not available AND local_path not found: {local_path}" + ) from e + + # Load local df + if p.suffix == ".pkl": + df = pd.read_pickle(p) + elif p.suffix == ".parquet": + df = pd.read_parquet(p) + else: + raise ValueError(f"Unsupported local dataset format: {p.suffix} (use .pkl or .parquet)") + + return _pandas_to_datasetdict(df, split="train") \ No newline at end of file diff --git a/engibench/problems/wings3D/templates/airfoil_analysis.py b/engibench/problems/wings3D/templates/wings3D_analysis.py similarity index 100% rename from engibench/problems/wings3D/templates/airfoil_analysis.py rename to engibench/problems/wings3D/templates/wings3D_analysis.py diff --git a/engibench/problems/wings3D/templates/airfoil_opt.py b/engibench/problems/wings3D/templates/wings3D_opt.py similarity index 100% rename from engibench/problems/wings3D/templates/airfoil_opt.py rename to engibench/problems/wings3D/templates/wings3D_opt.py diff --git a/engibench/problems/wings3D/utils.py b/engibench/problems/wings3D/utils.py index a1c6d2c6..6c1e6d11 100644 --- a/engibench/problems/wings3D/utils.py +++ b/engibench/problems/wings3D/utils.py @@ -1,10 +1,26 @@ -"""Utility functions for the airfoil problem.""" +"""Utility functions for the wings3D problem.""" import numpy as np import numpy.typing as npt import pandas as pd +def get_slice_coords(coords: npt.NDArray, slice_num: int) -> npt.NDArray[np.float32]: + """Extract one slice curve from coords. + + Expected shapes: + - (n_slices, n_points, 2): full wing sections, return coords[slice_num] + - (n_points, 2): already a single slice, return as-is + """ + arr = np.asarray(coords, dtype=np.float32) + + if arr.ndim == 3: + return arr[int(slice_num)] + if arr.ndim == 2: + return arr + + raise ValueError(f"Unexpected coords shape: {arr.shape}") + def _extract_connectivities(df_slice: pd.DataFrame) -> tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]: """Extract node connectivities from the dataframe slice. @@ -279,61 +295,51 @@ def reorder_coords(df_slice: pd.DataFrame) -> npt.NDArray[np.float32]: coords_x_reordered, coords_y_reordered, indices_reordered ) - return np.array([coords_x_reordered, coords_y_reordered]) + return np.column_stack((coords_x_reordered, coords_y_reordered)).astype(np.float32) - -def scale_coords( +def _scale_single_slice( coords: npt.NDArray[np.float64], blunted: bool = False, # noqa: FBT001, FBT002 xcut: float = 0.99, - min_trailing_edge_indices: float = 6, + min_trailing_edge_indices: int = 6, ) -> tuple[npt.NDArray[np.float64], bool]: - """Scales the coordinates to fit in the design space. + """Scale a single slice with shape (n_points, 2).""" + arr = np.asarray(coords, dtype=np.float64).copy() - Args: - coords (np.ndarray): The coordinates to scale. - blunted (bool): If True, the coordinates are assumed to be blunted. - xcut (float): The x coordinate of the cut, if the coordinates are blunted. - min_trailing_edge_indices (int): The minimum number of trailing edge indices to remove. + if arr.ndim != 2 or arr.shape[1] != 2: + raise ValueError(f"Expected single-slice coords with shape (n_points, 2), got {arr.shape}") - Returns: - np.ndarray: The scaled coordinates. - """ # Test if the coordinates are blunted or not - if not (blunted) and is_blunted(coords): + if (not blunted) and is_blunted(arr): blunted = True print( - "The coordinates may be blunted. However, blunted was not set to True. Will set blunted to True and continue, but please check the coordinates." + "The coordinates may be blunted. However, blunted was not set to True. " + "Will set blunted to True and continue, but please check the coordinates." ) - if not (blunted): - xcut = 1.0 - # Scale x coordinates to be xcut in length - airfoil_length = np.abs(np.max(coords[0, :]) - np.min(coords[0, :])) + airfoil_length = np.abs(np.max(arr[:, 0]) - np.min(arr[:, 0])) # Center the coordinates around the leading edge and scale them - coords[0, :] = xcut * (coords[0, :] - np.min(coords[0, :])) / airfoil_length - airfoil_length = np.abs(np.max(coords[0, :]) - np.min(coords[0, :])) + arr[:, 0] = xcut * (arr[:, 0] - np.min(arr[:, 0])) / airfoil_length # Shift the coordinates to be centered at 0 at the leading edge - leading_id = np.argmin(coords[0, :]) - y_dist = coords[1, leading_id] - coords[1, :] += -y_dist + leading_id = np.argmin(arr[:, 0]) + y_dist = arr[leading_id, 1] + arr[:, 1] += -y_dist + # Ensure the first and last points are the same - coords[0, 0] = xcut - coords[0, -1] = xcut - coords[1, -1] = coords[1, 0] - # Set the leading edge location + arr[0, 0] = xcut + arr[-1, 0] = xcut + arr[-1, 1] = arr[0, 1] if blunted: - coords_x = coords[0, :] - # Get all of the trailing edge indices, i.e where consecutive x coordinates are the same + coords_x = arr[:, 0] + err = 1e-5 x_gt = np.max(coords_x) * 0.99 trailing_edge_indices_l = np.where(np.abs(coords_x - np.roll(coords_x, -1)) < err)[0] trailing_edge_indices_r = np.where(np.abs(coords_x - np.roll(coords_x, 1)) < err)[0] - # Include any indices that are in either list trailing_edge_indices = np.unique(np.concatenate((trailing_edge_indices_l, trailing_edge_indices_r))) trailing_edge_indices = trailing_edge_indices[coords_x[trailing_edge_indices] >= x_gt] @@ -342,17 +348,52 @@ def scale_coords( while len(trailing_edge_indices) < min_trailing_edge_indices: trailing_edge_indices_l = np.where(np.abs(coords_x - np.roll(coords_x, -1)) < err)[0] trailing_edge_indices_r = np.where(np.abs(coords_x - np.roll(coords_x, 1)) < err)[0] - # Include any indices that are in either list trailing_edge_indices = np.unique(np.concatenate((trailing_edge_indices_l, trailing_edge_indices_r))) trailing_edge_indices = trailing_edge_indices[coords_x[trailing_edge_indices] >= x_gt] err *= 1.5 if err > err_stop: break - # Remove the trailing edge indices from the coordinates - coords = np.delete(coords, trailing_edge_indices[1:-1], axis=1) + if len(trailing_edge_indices) > 2: + arr = np.delete(arr, trailing_edge_indices[1:-1], axis=0) - return coords, blunted + return arr, blunted + + +def scale_coords( + coords: npt.NDArray[np.float64], + blunted: bool = False, # noqa: FBT001, FBT002 + xcut: float = 0.99, + min_trailing_edge_indices: int = 6, +) -> tuple[npt.NDArray[np.float64], bool]: + """Scales coordinates for either one slice or a full wing.""" + arr = np.asarray(coords, dtype=np.float64) + + if arr.ndim == 2: + return _scale_single_slice( + arr, + blunted=blunted, + xcut=xcut, + min_trailing_edge_indices=min_trailing_edge_indices, + ) + + if arr.ndim == 3: + scaled_slices = [] + any_blunted = blunted + + for curve in arr: + scaled_curve, curve_blunted = _scale_single_slice( + curve, + blunted=blunted, + xcut=xcut, + min_trailing_edge_indices=min_trailing_edge_indices, + ) + scaled_slices.append(scaled_curve) + any_blunted = any_blunted or curve_blunted + + return np.stack(scaled_slices), any_blunted + + raise ValueError(f"Unexpected coords shape: {arr.shape}") def calc_off_wall_distance( # noqa: PLR0913 @@ -392,38 +433,47 @@ def calc_off_wall_distance( # noqa: PLR0913 def is_blunted(coords: npt.NDArray[np.float64], delta_x_tol: float = 1e-5) -> bool: - """Checks if the coordinates are blunted or not. + """Checks if a single slice is blunted. Args: - coords (np.ndarray): The coordinates to check. + coords (np.ndarray): Slice coordinates with shape (n_points, 2). delta_x_tol (float): The tolerance for the x coordinate difference. Returns: bool: True if the coordinates are blunted, False otherwise. """ - # Check if the coordinates going away from the tip have a small delta y - coords_x = coords[0, :] - # Get all of the trailing edge indices, i.e where consecutive x coordinates are the same + arr = np.asarray(coords, dtype=np.float64) + + if arr.ndim != 2 or arr.shape[1] != 2: + raise ValueError(f"Expected single-slice coords with shape (n_points, 2), got {arr.shape}") + + coords_x = arr[:, 0] + x_gt = np.max(coords_x) * 0.99 trailing_edge_indices_l = np.where(np.abs(coords_x - np.roll(coords_x, -1)) < delta_x_tol)[0] trailing_edge_indices_r = np.where(np.abs(coords_x - np.roll(coords_x, 1)) < delta_x_tol)[0] - # Include any indices that are in either list + trailing_edge_indices = np.unique(np.concatenate((trailing_edge_indices_l, trailing_edge_indices_r))) trailing_edge_indices = trailing_edge_indices[coords_x[trailing_edge_indices] >= x_gt] - # check if we have no trailing edge indices return not len(trailing_edge_indices) <= 1 def calc_area(coords: npt.NDArray[np.float32]) -> float: - """Calculates the area of the airfoil. - - Args: - coords (np.ndarray): The coordinates of the airfoil. - - Returns: - float: The area of the airfoil. - """ - return 0.5 * np.absolute( - np.sum(coords[0, :] * np.roll(coords[1, :], -1)) - np.sum(coords[1, :] * np.roll(coords[0, :], -1)) - ) + """Calculates the area of a slice, or the summed area of all slices in a wing.""" + arr = np.asarray(coords, dtype=np.float32) + + if arr.ndim == 2: + x = arr[:, 0] + y = arr[:, 1] + return float(0.5 * np.absolute(np.sum(x * np.roll(y, -1)) - np.sum(y * np.roll(x, -1)))) + + if arr.ndim == 3: + total_area = 0.0 + for curve in arr: + x = curve[:, 0] + y = curve[:, 1] + total_area += 0.5 * np.absolute(np.sum(x * np.roll(y, -1)) - np.sum(y * np.roll(x, -1))) + return float(total_area) + + raise ValueError(f"Unexpected coords shape: {arr.shape}") \ No newline at end of file diff --git a/engibench/problems/wings3D/utils_wing3d.py b/engibench/problems/wings3D/utils_wing3d.py deleted file mode 100644 index ee4836ec..00000000 --- a/engibench/problems/wings3D/utils_wing3d.py +++ /dev/null @@ -1,26 +0,0 @@ -""" -Utility functions for the Wings3D problem. -""" - -from __future__ import annotations - -import numpy as np -import numpy.typing as npt - - -def get_slice_coords(coords: npt.NDArray, slice_num: int) -> npt.NDArray[np.float32]: - """ - Extract one slice curve from coords. - - Expected shapes: - - (9, 192, 2): full wing sections, return coords[slice_num] -> (192, 2) - - (192, 2): already a single slice, return as-is - """ - arr = np.asarray(coords, dtype=np.float32) - - if arr.ndim == 3: - return arr[int(slice_num)] - if arr.ndim == 2: - return arr - - raise ValueError(f"Unexpected coords shape: {arr.shape}") \ No newline at end of file diff --git a/engibench/problems/wings3D/v0.py b/engibench/problems/wings3D/v0.py index 303270bd..297e4d4b 100644 --- a/engibench/problems/wings3D/v0.py +++ b/engibench/problems/wings3D/v0.py @@ -1,4 +1,4 @@ -"""Airfoil problem. +"""Wings3D problem. Filename convention is that folder paths do not end with /. For example, /path/to/folder is correct, but /path/to/folder/ is not. """ @@ -6,7 +6,10 @@ import dataclasses from dataclasses import dataclass from dataclasses import field +from fileinput import filename from importlib.util import find_spec +from engibench.problems.power_electronics.utils import config +from engibench.problems.wings3D.dataset_hf_wings3d import load_wings3d_dataset import json import os import shutil @@ -19,6 +22,7 @@ import numpy as np import numpy.typing as npt import pandas as pd +from scipy.interpolate import interp1d from engibench.constraint import bounded from engibench.constraint import constraint @@ -27,18 +31,18 @@ from engibench.core import ObjectiveDirection from engibench.core import OptiStep from engibench.core import Problem -from engibench.problems.airfoil.pyopt_history import History -from engibench.problems.airfoil.templates import cli_interface -from engibench.problems.airfoil.utils import calc_area -from engibench.problems.airfoil.utils import calc_off_wall_distance -from engibench.problems.airfoil.utils import reorder_coords -from engibench.problems.airfoil.utils import scale_coords +from engibench.problems.wings3D.pyopt_history import History +from engibench.problems.wings3D.templates import cli_interface +from engibench.problems.wings3D.utils import calc_area +from engibench.problems.wings3D.utils import calc_off_wall_distance +from engibench.problems.wings3D.utils import reorder_coords +from engibench.problems.wings3D.utils import scale_coords from engibench.utils import container from engibench.utils.files import clone_dir # Allow loading pyoptsparse histories even if pyoptsparse is not installed: if find_spec("pyoptsparse") is None: - from engibench.problems.airfoil import fake_pyoptsparse + from engibench.problems.wings3D import fake_pyoptsparse sys.modules["pyoptsparse"] = fake_pyoptsparse @@ -73,17 +77,23 @@ def self_intersect(curve: npt.NDArray[np.float64]) -> tuple[int, npt.NDArray[np. @constraint(categories=IMPL) def does_not_self_intersect(design: DesignType) -> None: - """Check if a curve has no self intersections.""" - intersection = self_intersect(design["coords"]) - assert intersection is None, ( - f"design: Curve does self intersect at segment {intersection[0]}: {intersection[1]} -- {intersection[2]}" - ) + """Check that no wing slice has self intersections.""" + coords = np.asarray(design["coords"]) + + if coords.ndim != 3: + raise ValueError(f"Expected coords with shape (n_slices, n_points, 2), got {coords.shape}") + for slice_idx, curve in enumerate(coords): + intersection = self_intersect(curve) + assert intersection is None, ( + f"design: Slice {slice_idx} self-intersects at segment {intersection[0]}: " + f"{intersection[1]} -- {intersection[2]}" + ) -class Airfoil(Problem[DesignType]): - r"""Airfoil 2D shape optimization problem. +class Wings3D(Problem[DesignType]): + r"""Wings3D 3D shape optimization problem. - This problem simulates the performance of an airfoil in a 2D environment. An airfoil is represented by a set of 192 points that define its shape. The performance is evaluated by the [MACH-Aero](https://mdolab-mach-aero.readthedocs-hosted.com/en/latest/) simulator that computes the lift and drag coefficients of the airfoil. + This problem simulates the performance of a wing in a 3D environment. A wing is represented by a set of 192 points that define its shape. The performance is evaluated by the [MACH-Aero](https://mdolab-mach-aero.readthedocs-hosted.com/en/latest/) simulator that computes the lift and drag coefficients of the wing. """ version = 0 @@ -92,16 +102,20 @@ class Airfoil(Problem[DesignType]): ("cl", ObjectiveDirection.MAXIMIZE), ) + N_SLICES = 9 + N_POINTS = 192 + design_space = spaces.Dict( { - "coords": spaces.Box(low=0.0, high=1.0, shape=(2, 192), dtype=np.float32), + "coords": spaces.Box(low=0.0, high=1.0, shape=(N_SLICES, N_POINTS, 2), dtype=np.float32), "angle_of_attack": spaces.Box(low=0.0, high=10.0, shape=(1,), dtype=np.float32), } ) design_constraints = (does_not_self_intersect,) - dataset_id = "IDEALLab/airfoil_v0" + dataset_id = "Cashen/optiwing3d_engibench" container_id = "mdolab/public:u22-gcc-ompi-stable" __local_study_dir: str + #self.dataset_id = dataset_id @dataclass class Conditions: @@ -116,8 +130,8 @@ class Conditions: ] = 1e6 """Reynolds number""" area_initial: float = float("NAN") - """actual initial airfoil area""" - area_ratio_min: Annotated[float, bounded(lower=0.0, upper=1.2).category(THEORY)] = 0.7 + """actual initial wing area""" + area_ratio_min: float = 0.7 """Minimum ratio the initial area is allowed to decrease to i.e minimum_area = area_initial*area_target""" cl_target: float = 0.5 """Target lift coefficient to satisfy equality constraint""" @@ -154,7 +168,7 @@ def area_ratio_bound(area_ratio_min: float, area_initial: float, area_input_desi ) def __init__(self, seed: int = 0, base_directory: str | None = None) -> None: - """Initializes the Airfoil problem. + """Initializes the Wings3D problem. Args: seed (int): The random seed for the problem. @@ -166,7 +180,7 @@ def __init__(self, seed: int = 0, base_directory: str | None = None) -> None: self.__local_base_directory = base_directory else: self.__local_base_directory = os.getcwd() - self.__local_target_dir = self.__local_base_directory + "/engibench_studies/problems/airfoil" + self.__local_target_dir = self.__local_base_directory + "/engibench_studies/problems/wings3D" self.__local_template_dir = ( os.path.dirname(os.path.abspath(__file__)) + "/templates" ) # These templates are shipped with the lib @@ -175,7 +189,7 @@ def __init__(self, seed: int = 0, base_directory: str | None = None) -> None: # Docker target directory # This is used for files that are mounted into the docker container self.__docker_base_dir = "/home/mdolabuser/mount/engibench" - self.__docker_target_dir = self.__docker_base_dir + "/engibench_studies/problems/airfoil" + self.__docker_target_dir = self.__docker_base_dir + "/engibench_studies/problems/wings3D" super().__init__(seed=seed) @@ -186,7 +200,7 @@ def reset(self, seed: int | None = None, *, cleanup: bool = False) -> None: seed (int, optional): The seed to reset to. If None, a random seed is used. cleanup (bool): Deletes the previous study directory if True. """ - if cleanup: + if cleanup and os.path.exists(self.__local_study_dir): shutil.rmtree(self.__local_study_dir) super().reset(seed) @@ -243,8 +257,10 @@ def __design_to_simulator_input( ) # Save the design to a temporary file. Format to 1e-6 rounding - np.savetxt(self.__local_study_dir + "/" + filename + ".dat", scaled_design.transpose()) - + np.savetxt( + self.__local_study_dir + "/" + filename + ".dat", + scaled_design.reshape(-1, 2) + ) # Launches a docker container with the pre_process.py script # The script generates the mesh and FFD files bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && python {self.__docker_study_dir}/pre_process.py '{json.dumps(dataclasses.asdict(args))}'" @@ -306,59 +322,24 @@ def simulator_output_to_design(self, simulator_output: str | None = None) -> npt return reorder_coords(slice_df) - def simulate(self, design: DesignType, config: dict[str, Any] | None = None, mpicores: int = 4) -> npt.NDArray: - """Simulates the performance of an airfoil design. + def simulate(self, design: DesignType, config: dict[str, Any] | None = None, mpicores: int = 4): + """Simulate the performance (dummy for now).""" - Args: - design (dict): The design to simulate. - config (dict): A dictionary with configuration (e.g., boundary conditions, filenames) for the simulation. - mpicores (int): The number of MPI cores to use in the simulation. - - Returns: - dict: The performance of the design - each entry of the dict corresponds to a named objective value. - """ + # Ensure angle_of_attack is a float if isinstance(design["angle_of_attack"], np.ndarray): design["angle_of_attack"] = design["angle_of_attack"][0] - # pre-process the design and run the simulation + # Dummy simulation + drag = float(np.random.uniform(0.01, 0.05)) + lift = float(np.random.uniform(0.3, 1.2)) - # Prepares the airfoil_analysis.py script with the simulation configuration - conditions = self.Conditions() - config = config or {} - args = cli_interface.AnalysisParameters( - alpha=design["angle_of_attack"], - altitude=config.get("altitude", 10000), - temperature=config.get("temperature", 300), - reynolds=config.get("reynolds", conditions.reynolds), - mach=config.get("mach", conditions.mach), - use_altitude=config.get("use_altitude", False), - output_dir=config.get("output_dir", self.__docker_study_dir + "/output/"), - mesh_fname=config.get("mesh_fname", self.__docker_study_dir + "/design.cgns"), - task=cli_interface.Task[config["task"]] if "task" in config else cli_interface.Task.ANALYSIS, - ) - self.__design_to_simulator_input(design, mach=args.mach, reynolds=args.reynolds, temperature=args.temperature) - - # Launches a docker container with the airfoil_analysis.py script - # The script takes a mesh and ffd and performs an optimization - bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && mpirun -np {mpicores} python -m mpi4py {self.__docker_study_dir}/airfoil_analysis.py '{json.dumps(args.to_dict())}'" - assert self.container_id is not None, "Container ID is not set" - container.run( - command=["/bin/bash", "-c", bash_command], - image=self.container_id, - name="machaero", - mounts=[(self.__local_base_directory, self.__docker_base_dir)], - sync_uid=True, - ) - - outputs = np.load(self.__local_study_dir + "/output/outputs.npy") - lift = float(outputs[3]) - drag = float(outputs[4]) return np.array([drag, lift]) + def optimize( self, starting_point: DesignType, config: dict[str, Any] | None = None, mpicores: int = 4 ) -> tuple[DesignType, list[OptiStep]]: - """Optimizes the design of an airfoil. + """Optimizes the design of a wing. Args: starting_point (dict): The starting point for the optimization. @@ -374,7 +355,7 @@ def optimize( # pre-process the design and run the simulation filename = "candidate_design" - # Prepares the optimize_airfoil.py script with the optimization configuration + # Prepares the optimize_wings3D.py script with the optimization configuration fields = {f.name for f in dataclasses.fields(cli_interface.OptimizeParameters)} config = {key: val for key, val in (config or {}).items() if key in fields} if "area_initial" not in config: @@ -397,45 +378,16 @@ def optimize( **config, }, ) - self.__design_to_simulator_input( - starting_point, reynolds=args.reynolds, mach=args.reynolds, temperature=args.temperature, filename=filename - ) - - # Launches a docker container with the optimize_airfoil.py script - # The script takes a mesh and ffd and performs an optimization - bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && mpirun -np {mpicores} python -m mpi4py {self.__docker_study_dir}/airfoil_opt.py '{json.dumps(args.to_dict())}'" - assert self.container_id is not None, "Container ID is not set" - container.run( - command=["/bin/bash", "-c", bash_command], - image=self.container_id, - name="machaero", - mounts=[(self.__local_base_directory, self.__docker_base_dir)], - sync_uid=True, - ) + # Dummy optimization for now (skip simulation/CFD) + opt_design = { + "coords": starting_point["coords"], + "angle_of_attack": starting_point["angle_of_attack"], + } # post process -- extract the shape and objective values optisteps_history = [] - history = History(self.__local_study_dir + "/output/opt.hst") - call_counters = history.getCallCounters() - iters = list(map(int, call_counters)) if call_counters is not None else [] - - for i in range(len(iters)): - vals = history.read(int(iters[i])) - if vals is not None and "funcs" in vals and "obj" in vals["funcs"] and not vals["fail"]: - values = history.getValues(names=["obj"], callCounters=[i], allowSens=False, major=False, scale=True) - if values is not None and "obj" in values: - objective = values["obj"] - # flatten objective if it is a list - obj_np = np.array(objective) - if obj_np.ndim > 1: - obj_np = obj_np.flatten() - optisteps_history.append(OptiStep(obj_values=obj_np, step=vals["iter"])) - history.close() - - opt_coords = self.simulator_output_to_design() - - return {"coords": opt_coords, "angle_of_attack": starting_point["angle_of_attack"]}, optisteps_history + return opt_design, optisteps_history def render(self, design: DesignType, *, open_window: bool = False, save: bool = False) -> Figure: """Renders the design in a human-readable format. @@ -451,7 +403,11 @@ def render(self, design: DesignType, *, open_window: bool = False, save: bool = fig, ax = plt.subplots() coords = design["coords"] alpha = design["angle_of_attack"] - ax.scatter(coords[0], coords[1], s=10, alpha=0.7) + alpha_val = float(np.asarray(alpha).squeeze()) + + for curve in coords: + ax.plot(curve[:, 0], curve[:, 1], alpha=0.8) + ax.set_title(r"$\alpha$=" + str(np.round(alpha, 2)) + r"$^\circ$") ax.axis("equal") ax.axis("off") @@ -460,7 +416,7 @@ def render(self, design: DesignType, *, open_window: bool = False, save: bool = if open_window: plt.show() if save: - plt.savefig(self.__local_study_dir + "/airfoil.png", dpi=300, bbox_inches="tight") + plt.savefig(self.__local_study_dir + "/wings3D.png", dpi=300, bbox_inches="tight") plt.close(fig) return fig @@ -489,41 +445,76 @@ def render_optisteps(self, optisteps_history: list[OptiStep], *, open_window: bo plt.close(fig) return fig, ax - def random_design(self, dataset_split: str = "train", design_key: str = "initial_design") -> tuple[dict[str, Any], int]: - """Samples a valid random initial design. - - Args: - dataset_split (str): The key to use for the dataset. Defaults to "train". - design_key (str): The key to use for the design in the dataset. - Defaults to "initial_design". - - Returns: - tuple[dict[str, Any], int]: The valid random design and the index of the design in the dataset. - """ - rnd = self.np_random.integers(low=0, high=len(self.dataset[dataset_split][design_key]), dtype=int) - initial_design = self.dataset[dataset_split][design_key][rnd] - return {"coords": np.array(initial_design["coords"]), "angle_of_attack": initial_design["angle_of_attack"]}, rnd - + def random_design(self, dataset_split: str = "train") -> tuple[dict[str, Any], int]: + """Samples a valid random initial design.""" + dataset_split_data = self.dataset[dataset_split] + + rnd = 1 # self.np_random.integers(low=0, high=len(dataset_split_data["coords"]), dtype=int) + + raw_coords = dataset_split_data["coords"][rnd] + coords = np.stack([np.stack(s) for s in raw_coords]).astype(np.float32) + + # Remove duplicates and interpolate to ensure exactly N_POINTS per slice + from scipy.interpolate import interp1d + for i in range(len(coords)): + slice_coords = coords[i] + unique_coords, unique_indices = np.unique(slice_coords, axis=0, return_index=True) + unique_indices = np.sort(unique_indices) + unique_coords = slice_coords[unique_indices] + + if len(unique_coords) < self.N_POINTS: + # Interpolate to N_POINTS + t_old = np.linspace(0, 1, len(unique_coords)) + t_new = np.linspace(0, 1, self.N_POINTS) + interp_x = interp1d(t_old, unique_coords[:, 0], kind='linear') + interp_y = interp1d(t_old, unique_coords[:, 1], kind='linear') + new_x = interp_x(t_new) + new_y = interp_y(t_new) + slice_coords = np.column_stack((new_x, new_y)).astype(np.float32) + elif len(unique_coords) > self.N_POINTS: + slice_coords = unique_coords[:self.N_POINTS] + else: + slice_coords = unique_coords + + # # Check orientation and reverse if clockwise + # x = slice_coords[:, 0] + # y = slice_coords[:, 1] + # area = 0.5 * np.sum(x[:-1]*y[1:] - x[1:]*y[:-1]) + 0.5*(x[-1]*y[0] - x[0]*y[-1]) + # if area < 0: + # slice_coords = slice_coords[::-1] + + coords[i] = slice_coords + + angle_of_attack = dataset_split_data["alpha"][rnd] + + return { + "coords": coords, + "angle_of_attack": np.array([angle_of_attack], dtype=np.float32), + }, rnd if __name__ == "__main__": # Initialize the problem - problem = Airfoil(seed=0) + problem = Wings3D(seed=0) + problem.reset(seed=0) # Retrieve the dataset dataset = problem.dataset # Get random initial design and optimized conditions from the dataset + the index design, idx = problem.random_design() + print("Initial design (shape): ", design["coords"].shape) + print("Initial angle of attack: ", design["angle_of_attack"]) + + # Get the config conditions from the dataset + config = dataset["train"][idx] - # Get the config conditions from the dataset - config = dataset["train"].select_columns(problem.conditions_keys)[idx] # Simulate the design print("Simulation results: ", problem.simulate(design, config=config, mpicores=8)) # Cleanup the study directory; will delete the previous contents from simulate in this case - problem.reset(seed=1, cleanup=True) + #problem.reset(seed=1, cleanup=True) # Get design and conditions from the dataset, render design opt_design, optisteps_history = problem.optimize(design, config=config, mpicores=8) @@ -531,4 +522,4 @@ def random_design(self, dataset_split: str = "train", design_key: str = "initial print("Optimization history: ", optisteps_history) # Render the final optimized design - problem.render(opt_design, open_window=False, save=True) + #problem.render(opt_design, open_window=False, save=False) diff --git a/engibench/problems/wings3D/v0_offline.py b/engibench/problems/wings3D/v0_offline.py deleted file mode 100644 index e6964dec..00000000 --- a/engibench/problems/wings3D/v0_offline.py +++ /dev/null @@ -1,74 +0,0 @@ -""" -Wings3D (offline) dataset-backed problem. - -One sample corresponds to one slice (selected via slice_num) of a wing. -simulate() returns cd/cl stored in the dataset row. -""" - -from __future__ import annotations - -from typing import Any - -import numpy as np -import numpy.typing as npt -from gymnasium import spaces -from matplotlib.figure import Figure -import matplotlib.pyplot as plt - -from engibench.core import ObjectiveDirection, Problem -from engibench.problems.wings3D.dataset_hf_wings3d import load_wings3d_dataset -from engibench.problems.wings3D.utils_wing3d import get_slice_coords - -DesignType = dict[str, Any] - - -class Wings3DOffline(Problem[DesignType]): - version = 0 - objectives = (("cd", ObjectiveDirection.MINIMIZE), ("cl", ObjectiveDirection.MAXIMIZE)) - - design_space = spaces.Dict( - { - "coords": spaces.Box(low=-1e6, high=1e6, shape=(192, 2), dtype=np.float32), - "alpha": spaces.Box(low=-90.0, high=90.0, shape=(1,), dtype=np.float32), - } - ) - - dataset_id: str = "IDEALLab/wings3d_v0" # replace with real HF id when available - - def __init__(self, seed: int = 0) -> None: - super().__init__(seed=seed) - self.dataset = load_wings3d_dataset(self.dataset_id) - - def random_design(self, dataset_split: str = "train") -> tuple[DesignType, dict[str, Any], int]: - ds = self.dataset[dataset_split] - idx = int(self.np_random.integers(low=0, high=len(ds), dtype=int)) - row = ds[idx] - - slice_num = int(row["slice_num"]) - coords_slice = get_slice_coords(row["coords"], slice_num) - alpha = np.asarray([row.get("alpha", 0.0)], dtype=np.float32) - - design: DesignType = {"coords": coords_slice, "alpha": alpha} - config = {"split": dataset_split, "idx": idx} - return design, config, idx - - def simulate(self, design: DesignType, config: dict[str, Any] | None = None, **kwargs) -> npt.NDArray[np.float64]: - if not config or "split" not in config or "idx" not in config: - raise ValueError("simulate() requires config with keys {'split','idx'} (use random_design()).") - - row = self.dataset[config["split"]][int(config["idx"])] - cd = float(row.get("cd_val", np.nan)) - cl = float(row.get("cl_val", np.nan)) - return np.array([cd, cl], dtype=np.float64) - - def render(self, design: DesignType, *, open_window: bool = False, save: bool = False) -> Figure: - coords = np.asarray(design["coords"], dtype=np.float32) - fig, ax = plt.subplots() - ax.plot(coords[:, 0], coords[:, 1]) - ax.axis("equal") - ax.axis("off") - if open_window: - plt.show() - plt.close(fig) - return fig - \ No newline at end of file diff --git a/engibench/utils/all_problems.py b/engibench/utils/all_problems.py index 87ecb9a7..93359bcf 100644 --- a/engibench/utils/all_problems.py +++ b/engibench/utils/all_problems.py @@ -10,26 +10,45 @@ def list_problems(base_module: Any = engibench.problems) -> dict[str, type[Problem]]: - """Return a dict containing all available Problem instances defined in submodules of `base_module`.""" + """Return a dict containing all available Problem classes defined in submodules of `base_module`.""" module_path = os.path.dirname(base_module.__file__) modules = pkgutil.iter_modules(path=[module_path], prefix="") - return {m.name: extract_problem(importlib.import_module(f"{base_module.__package__}.{m.name}")) for m in modules} + out: dict[str, type[Problem]] = {} + for m in modules: + modname = f"{base_module.__package__}.{m.name}" + try: + module = importlib.import_module(modname) + except (ModuleNotFoundError, ImportError): + # Optional dependency missing (e.g., ceviche) or other import error + continue -def extract_problem(module: Any) -> type[Problem]: + p = extract_problem(module) + if p is None: + continue + + out[m.name] = p + + return out + + +def extract_problem(module: Any) -> type[Problem] | None: """Get a `Problem` class defined in a module. + Returns None if the module contains no `Problem` classes. Raises an exception if the module contains multiple `Problem` classes. """ problem_types = [ o for o in vars(module).values() if isinstance(o, type) and issubclass(o, Problem) and o is not Problem ] - try: - (p,) = problem_types - except ValueError: + + if len(problem_types) == 0: + return None + + if len(problem_types) > 1: msg = f"Only one problem per module is allowed. Got {', '.join(p.__name__ for p in problem_types)}" - raise ValueError(msg) from None - return p + raise ValueError(msg) + return problem_types[0] BUILTIN_PROBLEMS = list_problems() From 5ae47f97549805a5a934211d1e6a165f9e01fb26 Mon Sep 17 00:00:00 2001 From: Anna Delbeke Date: Mon, 16 Mar 2026 15:13:32 +0100 Subject: [PATCH 09/10] Fix ruff and mypy issues --- engibench/problems/wings3D/v0.py | 66 ++++++++++++-------------------- 1 file changed, 25 insertions(+), 41 deletions(-) diff --git a/engibench/problems/wings3D/v0.py b/engibench/problems/wings3D/v0.py index 297e4d4b..5d3b5887 100644 --- a/engibench/problems/wings3D/v0.py +++ b/engibench/problems/wings3D/v0.py @@ -6,10 +6,7 @@ import dataclasses from dataclasses import dataclass from dataclasses import field -from fileinput import filename from importlib.util import find_spec -from engibench.problems.power_electronics.utils import config -from engibench.problems.wings3D.dataset_hf_wings3d import load_wings3d_dataset import json import os import shutil @@ -31,7 +28,6 @@ from engibench.core import ObjectiveDirection from engibench.core import OptiStep from engibench.core import Problem -from engibench.problems.wings3D.pyopt_history import History from engibench.problems.wings3D.templates import cli_interface from engibench.problems.wings3D.utils import calc_area from engibench.problems.wings3D.utils import calc_off_wall_distance @@ -48,6 +44,9 @@ DesignType = dict[str, Any] +# Expected dimensionality of the wing coordinate array: (n_slices, n_points, 2) +COORDS_NDIMS = 3 + def self_intersect(curve: npt.NDArray[np.float64]) -> tuple[int, npt.NDArray[np.float64], npt.NDArray[np.float64]] | None: """Determines if two segments a and b intersect.""" @@ -80,7 +79,7 @@ def does_not_self_intersect(design: DesignType) -> None: """Check that no wing slice has self intersections.""" coords = np.asarray(design["coords"]) - if coords.ndim != 3: + if coords.ndim != COORDS_NDIMS: raise ValueError(f"Expected coords with shape (n_slices, n_points, 2), got {coords.shape}") for slice_idx, curve in enumerate(coords): @@ -90,6 +89,7 @@ def does_not_self_intersect(design: DesignType) -> None: f"{intersection[1]} -- {intersection[2]}" ) + class Wings3D(Problem[DesignType]): r"""Wings3D 3D shape optimization problem. @@ -115,7 +115,6 @@ class Wings3D(Problem[DesignType]): dataset_id = "Cashen/optiwing3d_engibench" container_id = "mdolab/public:u22-gcc-ompi-stable" __local_study_dir: str - #self.dataset_id = dataset_id @dataclass class Conditions: @@ -257,10 +256,7 @@ def __design_to_simulator_input( ) # Save the design to a temporary file. Format to 1e-6 rounding - np.savetxt( - self.__local_study_dir + "/" + filename + ".dat", - scaled_design.reshape(-1, 2) - ) + np.savetxt(self.__local_study_dir + "/" + filename + ".dat", scaled_design.reshape(-1, 2)) # Launches a docker container with the pre_process.py script # The script generates the mesh and FFD files bash_command = f"source /home/mdolabuser/.bashrc_mdolab && cd {self.__docker_base_dir} && python {self.__docker_study_dir}/pre_process.py '{json.dumps(dataclasses.asdict(args))}'" @@ -324,18 +320,19 @@ def simulator_output_to_design(self, simulator_output: str | None = None) -> npt def simulate(self, design: DesignType, config: dict[str, Any] | None = None, mpicores: int = 4): """Simulate the performance (dummy for now).""" + del config, mpicores # Ensure angle_of_attack is a float if isinstance(design["angle_of_attack"], np.ndarray): design["angle_of_attack"] = design["angle_of_attack"][0] # Dummy simulation - drag = float(np.random.uniform(0.01, 0.05)) - lift = float(np.random.uniform(0.3, 1.2)) + rng = np.random.default_rng() + drag = float(rng.uniform(0.01, 0.05)) + lift = float(rng.uniform(0.3, 1.2)) return np.array([drag, lift]) - def optimize( self, starting_point: DesignType, config: dict[str, Any] | None = None, mpicores: int = 4 ) -> tuple[DesignType, list[OptiStep]]: @@ -349,6 +346,8 @@ def optimize( Returns: tuple[dict[str, Any], list[OptiStep]]: The optimized design and its performance. """ + del mpicores + if isinstance(starting_point["angle_of_attack"], np.ndarray): starting_point["angle_of_attack"] = starting_point["angle_of_attack"][0] @@ -362,7 +361,7 @@ def optimize( raise ValueError("optimize(): config is missing the required parameter 'area_initial'") if "opt" in config: config["opt"] = cli_interface.Algorithm[config["opt"]] - args = cli_interface.OptimizeParameters( + _args = cli_interface.OptimizeParameters( **{ **dataclasses.asdict(self.Conditions()), "alpha": starting_point["angle_of_attack"], @@ -385,7 +384,7 @@ def optimize( } # post process -- extract the shape and objective values - optisteps_history = [] + optisteps_history: list[OptiStep] = [] return opt_design, optisteps_history @@ -402,8 +401,7 @@ def render(self, design: DesignType, *, open_window: bool = False, save: bool = """ fig, ax = plt.subplots() coords = design["coords"] - alpha = design["angle_of_attack"] - alpha_val = float(np.asarray(alpha).squeeze()) + alpha = float(np.asarray(design["angle_of_attack"]).squeeze()) for curve in coords: ax.plot(curve[:, 0], curve[:, 1], alpha=0.8) @@ -455,35 +453,26 @@ def random_design(self, dataset_split: str = "train") -> tuple[dict[str, Any], i coords = np.stack([np.stack(s) for s in raw_coords]).astype(np.float32) # Remove duplicates and interpolate to ensure exactly N_POINTS per slice - from scipy.interpolate import interp1d + for i in range(len(coords)): slice_coords = coords[i] unique_coords, unique_indices = np.unique(slice_coords, axis=0, return_index=True) unique_indices = np.sort(unique_indices) unique_coords = slice_coords[unique_indices] - + if len(unique_coords) < self.N_POINTS: # Interpolate to N_POINTS t_old = np.linspace(0, 1, len(unique_coords)) t_new = np.linspace(0, 1, self.N_POINTS) - interp_x = interp1d(t_old, unique_coords[:, 0], kind='linear') - interp_y = interp1d(t_old, unique_coords[:, 1], kind='linear') + interp_x = interp1d(t_old, unique_coords[:, 0], kind="linear") + interp_y = interp1d(t_old, unique_coords[:, 1], kind="linear") new_x = interp_x(t_new) new_y = interp_y(t_new) slice_coords = np.column_stack((new_x, new_y)).astype(np.float32) elif len(unique_coords) > self.N_POINTS: - slice_coords = unique_coords[:self.N_POINTS] + slice_coords = unique_coords[: self.N_POINTS] else: slice_coords = unique_coords - - # # Check orientation and reverse if clockwise - # x = slice_coords[:, 0] - # y = slice_coords[:, 1] - # area = 0.5 * np.sum(x[:-1]*y[1:] - x[1:]*y[:-1]) + 0.5*(x[-1]*y[0] - x[0]*y[-1]) - # if area < 0: - # slice_coords = slice_coords[::-1] - - coords[i] = slice_coords angle_of_attack = dataset_split_data["alpha"][rnd] @@ -492,6 +481,7 @@ def random_design(self, dataset_split: str = "train") -> tuple[dict[str, Any], i "angle_of_attack": np.array([angle_of_attack], dtype=np.float32), }, rnd + if __name__ == "__main__": # Initialize the problem @@ -506,20 +496,14 @@ def random_design(self, dataset_split: str = "train") -> tuple[dict[str, Any], i print("Initial design (shape): ", design["coords"].shape) print("Initial angle of attack: ", design["angle_of_attack"]) - # Get the config conditions from the dataset - config = dataset["train"][idx] - + # Get the config conditions from the dataset + dataset_config = dataset["train"][idx] # Simulate the design - print("Simulation results: ", problem.simulate(design, config=config, mpicores=8)) - - # Cleanup the study directory; will delete the previous contents from simulate in this case - #problem.reset(seed=1, cleanup=True) + print("Simulation results: ", problem.simulate(design, config=dataset_config, mpicores=8)) # Get design and conditions from the dataset, render design - opt_design, optisteps_history = problem.optimize(design, config=config, mpicores=8) + opt_design, optisteps_history = problem.optimize(design, config=dataset_config, mpicores=8) print("Optimized design: ", opt_design) print("Optimization history: ", optisteps_history) - # Render the final optimized design - #problem.render(opt_design, open_window=False, save=False) From 3e209844bdd1048aa91aaf758ea1112cdaf4e600 Mon Sep 17 00:00:00 2001 From: Anna Delbeke Date: Tue, 17 Mar 2026 15:27:48 +0100 Subject: [PATCH 10/10] add wings3D problem --- engibench/problems/wings3D/v0.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/engibench/problems/wings3D/v0.py b/engibench/problems/wings3D/v0.py index 5d3b5887..d57abbbc 100644 --- a/engibench/problems/wings3D/v0.py +++ b/engibench/problems/wings3D/v0.py @@ -474,7 +474,7 @@ def random_design(self, dataset_split: str = "train") -> tuple[dict[str, Any], i else: slice_coords = unique_coords - angle_of_attack = dataset_split_data["alpha"][rnd] + angle_of_attack = dataset_split_data["alpha"][rnd]euler-tunnel start return { "coords": coords,