diff --git a/engibench/problems/__init__.py b/engibench/problems/__init__.py index 069a7ff8..c619c079 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 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/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_hf_wings3d.py b/engibench/problems/wings3D/dataset_hf_wings3d.py new file mode 100644 index 00000000..cbcb3d6e --- /dev/null +++ b/engibench/problems/wings3D/dataset_hf_wings3d.py @@ -0,0 +1,54 @@ +""" +Dataset loader for the Wings3D problem. + +Tries Hugging Face first. Optionally falls back to a local file for development. +""" + +from __future__ import annotations + +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, 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: + 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/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/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/templates/wings3D_analysis.py b/engibench/problems/wings3D/templates/wings3D_analysis.py new file mode 100644 index 00000000..e73b1ebd --- /dev/null +++ b/engibench/problems/wings3D/templates/wings3D_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/wings3D_opt.py b/engibench/problems/wings3D/templates/wings3D_opt.py new file mode 100644 index 00000000..6a0b9c75 --- /dev/null +++ b/engibench/problems/wings3D/templates/wings3D_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/utils.py b/engibench/problems/wings3D/utils.py new file mode 100644 index 00000000..6c1e6d11 --- /dev/null +++ b/engibench/problems/wings3D/utils.py @@ -0,0 +1,479 @@ +"""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. + + 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.column_stack((coords_x_reordered, coords_y_reordered)).astype(np.float32) + +def _scale_single_slice( + 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]: + """Scale a single slice with shape (n_points, 2).""" + arr = np.asarray(coords, dtype=np.float64).copy() + + if arr.ndim != 2 or arr.shape[1] != 2: + raise ValueError(f"Expected single-slice coords with shape (n_points, 2), got {arr.shape}") + + # Test if the coordinates are blunted or not + 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." + ) + + # Scale x coordinates to be xcut in length + airfoil_length = np.abs(np.max(arr[:, 0]) - np.min(arr[:, 0])) + + # Center the coordinates around the leading edge and scale them + 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(arr[:, 0]) + y_dist = arr[leading_id, 1] + arr[:, 1] += -y_dist + + # Ensure the first and last points are the same + arr[0, 0] = xcut + arr[-1, 0] = xcut + arr[-1, 1] = arr[0, 1] + + if blunted: + 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] + 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] + 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 + + if len(trailing_edge_indices) > 2: + arr = np.delete(arr, trailing_edge_indices[1:-1], axis=0) + + 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 + 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 a single slice is blunted. + + Args: + 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. + """ + 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] + + 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] + + return not len(trailing_edge_indices) <= 1 + + +def calc_area(coords: npt.NDArray[np.float32]) -> float: + """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/v0.py b/engibench/problems/wings3D/v0.py new file mode 100644 index 00000000..d57abbbc --- /dev/null +++ b/engibench/problems/wings3D/v0.py @@ -0,0 +1,509 @@ +"""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. +""" + +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 scipy.interpolate import interp1d + +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.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.wings3D import fake_pyoptsparse + + sys.modules["pyoptsparse"] = fake_pyoptsparse + +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.""" + # 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 that no wing slice has self intersections.""" + coords = np.asarray(design["coords"]) + + 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): + 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 Wings3D(Problem[DesignType]): + r"""Wings3D 3D shape optimization problem. + + 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 + objectives: tuple[tuple[str, ObjectiveDirection], ...] = ( + ("cd", ObjectiveDirection.MINIMIZE), + ("cl", ObjectiveDirection.MAXIMIZE), + ) + + N_SLICES = 9 + N_POINTS = 192 + + design_space = spaces.Dict( + { + "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 = "Cashen/optiwing3d_engibench" + 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 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""" + + 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 Wings3D 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/wings3D" + 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/wings3D" + + 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 and os.path.exists(self.__local_study_dir): + 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.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))}'" + 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): + """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 + 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]]: + """Optimizes the design of a wing. + + 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. + """ + del mpicores + + 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_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: + 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, + }, + ) + # 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: list[OptiStep] = [] + + 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. + + 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 = float(np.asarray(design["angle_of_attack"]).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") + ax.set_xlim((-0.005, 1.005)) + + if open_window: + plt.show() + if save: + plt.savefig(self.__local_study_dir + "/wings3D.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") -> 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 + + 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 + + angle_of_attack = dataset_split_data["alpha"][rnd]euler-tunnel start + + return { + "coords": coords, + "angle_of_attack": np.array([angle_of_attack], dtype=np.float32), + }, rnd + + +if __name__ == "__main__": + # Initialize the problem + + 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 + dataset_config = dataset["train"][idx] + + # Simulate the design + 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=dataset_config, mpicores=8) + print("Optimized design: ", opt_design) + print("Optimization history: ", optisteps_history) + 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()