From 23b8c5dd8fab815e672d1b9edfdb107a35b2faf3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Thu, 8 Jan 2026 16:03:18 +0100 Subject: [PATCH 1/3] README.md: Add missing comma in bibliography --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index f0d24dc1..008136fd 100644 --- a/README.md +++ b/README.md @@ -111,7 +111,7 @@ If you use EngiBench in your research, please cite the following paper: url = {https://openreview.net/forum?id=YowD33Q89V}, urldate = {2025-10-07}, author = {Felten, Florian and Apaza, Gabriel and B\¨aunlich, Gerhard and Diniz, Cashen and Dong, Xuliang and Drake, Arthur and Habibi, Milad and Hoffman, Nathaniel J. and Keeler, Matthew and Massoudi, Soheyl and VanGessel, Francis G. and Fuge, Mark}, - booktitle = {Proceedings of the 39th Conference on Neural Information Processing Systems ({NeurIPS} 2025)} + booktitle = {Proceedings of the 39th Conference on Neural Information Processing Systems ({NeurIPS} 2025)}, year = {2025}, } ``` From 8e16f7105202b40decff8e601d885037b83bd991 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Mon, 26 Jan 2026 11:22:30 +0100 Subject: [PATCH 2/3] docs: Fix invalid refs and unincluded content --- docs/conf.py | 12 +++++++----- docs/index.md | 5 +---- docs/tutorials/new_problem.md | 8 ++++---- 3 files changed, 12 insertions(+), 13 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 5fc8ef99..bdefe203 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -10,11 +10,7 @@ # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. # -# import os -# import sys -# sys.path.insert(0, os.path.abspath('.')) -# -- Project information ----------------------------------------------------- import os from pathlib import Path import sys @@ -52,7 +48,7 @@ # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. # This pattern also affects html_static_path and html_extra_path. -exclude_patterns: list[str] = [] +exclude_patterns: list[str] = ["README.md"] # Napoleon settings napoleon_use_ivar = True @@ -131,3 +127,9 @@ os.path.realpath(os.path.join(os.path.dirname(__file__), "..", "tests", "tools")), os.path.realpath(os.path.join(os.path.dirname(__file__), "utils")), ] + +myst_url_schemes = { + "http": None, + "https": None, + "source": "https://github.com/IDEALLab/EngiBench/tree/main/{{path}}", +} diff --git a/docs/index.md b/docs/index.md index cf72bec4..4dbafdc2 100644 --- a/docs/index.md +++ b/docs/index.md @@ -69,11 +69,8 @@ EngiOpt ``` ```{toctree} -:hidden: -:caption: Utils -utils/container -utils/slurm +utils/index ``` ```{toctree} diff --git a/docs/tutorials/new_problem.md b/docs/tutorials/new_problem.md index 1ebb95dc..d048ec00 100644 --- a/docs/tutorials/new_problem.md +++ b/docs/tutorials/new_problem.md @@ -12,7 +12,7 @@ Note: ### Code In general, follow the `beams2d/` example. -1. Create a new problem module in [engibench/problems/](engibench/problems/) following the following layout (e.g. [engibench/problems/beams2d/](engibench/problems/beams2d/)), where you later also can add other versions / variant of the problem: +1. Create a new problem module in [engibench/problems/](source:engibench/problems/) following the following layout (e.g. [engibench/problems/beams2d/](source:engibench/problems/beams2d/)), where you later also can add other versions / variant of the problem: ``` 📦 engibench └─ 📂 problems @@ -35,7 +35,7 @@ In general, follow the `beams2d/` example. Ideally, all non-breaking changes should not create a new versioned module. Also in many cases, code duplication can be avoided, by introducing a new parameter to the problem class. -2. Define your problem class that implements the `Problem` interface with its functions and attributes in `problems/new_problem/v0.py` (e.g. [beams2d/v0.py](engibench/problems/beams2d/v0.py)). +2. Define your problem class that implements the `Problem` interface with its functions and attributes in `problems/new_problem/v0.py` (e.g. [beams2d/v0.py](source:engibench/problems/beams2d/v0.py)). `problems/new_problem/v0.py` ```py @@ -52,8 +52,8 @@ In general, follow the `beams2d/` example. #### Documentation 1. Install necessary documentation tools: `pip install ".[doc]"`. -2. If it is a new problem family, add a new `.md` file in [docs/problems/](docs/problems/) following - the existing structure and add your problem family in the `toctree` of [docs/problems/index.md](docs/problems/index.md). +2. If it is a new problem family, add a new `.md` file in [docs/problems/](source:docs/problems/) following + the existing structure and add your problem family in the `toctree` of [docs/problems/index.md](source:docs/problems/index.md). 3. Add a problem markdown file to the `toctree` in `docs/problems/new_problem.md`. In the md file, use EngiBench's own `problem` directive: ``````md # Your Problem From e26abeb643d484e5173855d5b5a26dd5bd65cdf5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Gerhard=20Br=C3=A4unlich?= Date: Thu, 8 Jan 2026 16:10:32 +0100 Subject: [PATCH 3/3] docs: Improve generation of problem docs --- docs/README.md | 30 ++- docs/_ext/problem_doc.py | 230 +++++++++++++++++++-- docs/problems/airfoil.md | 162 ++++++++++++++- docs/problems/beams2d.md | 103 ++++++++- docs/problems/heatconduction2d.md | 38 +++- docs/problems/heatconduction3d.md | 44 +++- docs/problems/photonics2d.md | 148 ++++++++++++- docs/problems/power_electronics.md | 79 ++++++- docs/problems/thermoelastic2d.md | 52 ++++- docs/tutorials/new_problem.md | 35 +++- engibench/problems/airfoil/v0.py | 168 +-------------- engibench/problems/beams2d/v0.py | 110 +--------- engibench/problems/beams2d/v1.py | 3 +- engibench/problems/heatconduction2d/v0.py | 42 +--- engibench/problems/heatconduction3d/v0.py | 49 +---- engibench/problems/photonics2d/v0.py | 155 +------------- engibench/problems/power_electronics/v0.py | 80 ------- engibench/problems/thermoelastic2d/v0.py | 66 +----- 18 files changed, 907 insertions(+), 687 deletions(-) diff --git a/docs/README.md b/docs/README.md index 6dee3181..3f84a516 100644 --- a/docs/README.md +++ b/docs/README.md @@ -11,16 +11,38 @@ Fork EngiBench and edit the docstring in the problem's Python file. Then, pip in ### Adding a new problem Ensure the problem is in EngiBench (or your fork). Ensure that its Python file has a properly formatted markdown docstring. Install using `pip install -e .[doc]` and add a markdown file in the [./problems/](problems/) of the repo. -Use Engibench's own `problem` directive in the docs of your new problem: +Use EngiBench's own `problem:table` and `problem:conditions` directives in the docs of your new problem: ``````md # Your Problem - ``` {problem} your_problem + ``` {problem:table} + :lead: Chuck Norris @chucknorris ``` + + ... + + ## Conditions + + ``` {problem:conditions} + ``` + + ... `````` -Here, `your_problem` must match the name of the module where your problem class is defined. -This will automatically include the docstrings of your `Problem` class as well as a table with its metadata. Then complete the [other steps](#other-steps). +**`problem:table`**: This directive extracts metadata from a problem and +inserts a table filled with the metadata. +By default, the directive will try to import the problem `engibench.problems.`, where `` is the filename (without `.md` extension) of the markdown file the directive is used. + +Options (optional): +* `:problem_id:` override ``, +* `:lead:` Add a row "Lead" to the table, containing the specified value. + If the value ends with `@username`, a link to `https://github.com/username` will be inserted. + +**`problem:conditions`**: This directive lists the conditions extracted from a problem as in the "Conditions" row produced by the `problem:table` directive. The `` is determined the same way as in `problem:table`. + +Options: +* `:problem_id:` override ``, +* `:defaults:` include default values in the list of conditions #### Other steps diff --git a/docs/_ext/problem_doc.py b/docs/_ext/problem_doc.py index c26cfe85..e94a5861 100644 --- a/docs/_ext/problem_doc.py +++ b/docs/_ext/problem_doc.py @@ -8,7 +8,7 @@ ````` """ -from collections.abc import Iterator, Sequence +from collections.abc import Iterable, Iterator, Sequence import contextlib import dataclasses import importlib.abc @@ -16,34 +16,66 @@ import inspect import sys from types import ModuleType -from typing import Any +from typing import Any, ClassVar, get_type_hints, TYPE_CHECKING import unittest.mock from docutils import nodes +from docutils.parsers.rst import directives +from sphinx import addnodes from sphinx.application import Sphinx +from sphinx.domains import Domain from sphinx.util.docutils import SphinxDirective +from sphinx.util.typing import ExtensionMetadata + +if TYPE_CHECKING: + from sphinx.builders import Builder + from sphinx.environment import BuildEnvironment MODULE_WHITELIST = frozenset(["engibench"]) MODULE_EXTRA_MEMBERS = {"networkx": ["Graph"], "gymnasium": ["spaces"]} -def setup(app: Sphinx) -> None: +def setup(app: Sphinx) -> ExtensionMetadata: """Add extension to sphinx.""" - app.add_directive("problem", ProblemDirective) + app.add_domain(ProblemDomain) + return { + "version": "0.1", + "parallel_read_safe": True, + "parallel_write_safe": True, + } + + +class Lead: + """Option to specify a lead in the problem directive.""" + + caption = "Lead" + def __init__(self, value: str) -> None: + self.handle = None + self.name, self.handle = value.split(" @", 1) if " @" in value else (value, None) -class ProblemDirective(SphinxDirective): - required_arguments = 1 + def to_node(self) -> nodes.Node: + p = nodes.Text(self.name) + if self.handle: + node = nodes.paragraph() + node += [ + p, + nodes.Text(" "), + nodes.reference(refuri=f"https://github.com/{self.handle}", text="@" + self.handle), + ] + return node + return p + + +class ProblemTableDirective(SphinxDirective): + option_spec: ClassVar[dict[str, Any]] = {"lead": Lead, "problem_id": str} def run(self) -> list[Any]: - with mock_imports(MODULE_WHITELIST, extra_members=MODULE_EXTRA_MEMBERS): - from engibench.core import ObjectiveDirection - from engibench.utils.all_problems import BUILTIN_PROBLEMS + problem_id = self.options.get("problem_id") or problem_id_from_docname(self.env.docname) + problem = import_problem(problem_id) + ObjectiveDirection = import_objective_direction() # noqa: N806 - problem_id = self.arguments[0].strip() - problem = BUILTIN_PROBLEMS[problem_id] - docstring = unindent(problem.__doc__) if problem.__doc__ is not None else "" - docstring = inspect.cleandoc(docstring) + docstring = inspect.getdoc(problem) image = nodes.image(uri=f"../_static/img/problems/{problem_id}.png", width="450px", align="center") @@ -51,16 +83,19 @@ def run(self) -> list[Any]: f"{obj}: ↑" if direction == ObjectiveDirection.MAXIMIZE else f"{obj}: ↓" for obj, direction in problem.objectives ] - conditions = [f"{f.name}: {f.default}" for f in dataclasses.fields(problem.Conditions)] + conditions = read_dataclass(problem.Conditions) + + lead = self.options.get("lead") tab_data = [ ("Version", str(problem.version)), ("Design space", make_code(repr(problem.design_space))), ("Objectives", make_multiline(objectives)), - ("Conditions", make_multiline(conditions)), + ("Conditions", make_simple_field_list(conditions)), ("Dataset", make_link(problem.dataset_id, f"https://huggingface.co/datasets/{problem.dataset_id}")), ("Container", make_code(problem.container_id) if problem.container_id is not None else None), ("Import", make_code(f"from {problem.__module__} import {problem.__name__}")), + *([("Lead", lead.to_node())] if lead is not None else []), ] # Very ugly hack to retain the order of children @@ -73,15 +108,60 @@ def run(self) -> list[Any]: sec.clear() sec.extend(header) - return [image, make_table(tab_data), *body] + return [image, *body, make_table(tab_data)] + + +def problem_id_from_docname(docname: str) -> str: + _, problem_id = docname.rsplit("/", 1) + return problem_id + + +class ConditionsDirective(SphinxDirective): + option_spec: ClassVar[dict[str, Any]] = {"problem_id": str, "defaults": bool} + + def run(self) -> list[Any]: + problem_id = self.options.get("problem_id") or problem_id_from_docname(self.env.docname) + problem = import_problem(problem_id) + + conditions = read_dataclass(problem.Conditions) + return [make_simple_field_list(conditions, defaults=self.options.get("defaults", False))] + + +class ProblemDomain(Domain): + name = "problem" + label = "Engibench Problem" + + directives: ClassVar[dict[str, SphinxDirective]] = { + "table": ProblemTableDirective, + "conditions": ConditionsDirective, + } + + def resolve_any_xref( # noqa: PLR0913 + self, + env: "BuildEnvironment", # noqa: ARG002 + fromdocname: str, # noqa: ARG002 + builder: "Builder", # noqa: ARG002 + target: str, # noqa: ARG002 + node: addnodes.pending_xref, # noqa: ARG002 + contnode: nodes.Element, # noqa: ARG002 + ) -> list[tuple[str, nodes.reference]]: + return [] + + +def import_objective_direction() -> type[Any]: + """Import the ObjectiveDirection enum without requiring engibench dependencies.""" + with mock_imports(MODULE_WHITELIST, extra_members=MODULE_EXTRA_MEMBERS): + from engibench.core import ObjectiveDirection # noqa: PLC0415 + return ObjectiveDirection -def make_section(title: str, section_id: str, body: list[Any]) -> nodes.section: - sec = nodes.section(ids=[section_id]) - sec += nodes.title(text=title) - for element in body: - sec += element - return sec + +def import_problem(problem_id: str) -> Any: + """Import problem metadata without requiring engibench dependencies.""" + with mock_imports(MODULE_WHITELIST, extra_members=MODULE_EXTRA_MEMBERS): + from engibench.utils.all_problems import BUILTIN_PROBLEMS # noqa: PLC0415 + + return BUILTIN_PROBLEMS[problem_id] def make_link(text: str, uri: str) -> nodes.paragraph: @@ -118,6 +198,59 @@ def make_table(tab_data: list[tuple[str, Any]]) -> nodes.table: return table +@dataclasses.dataclass +class Field: + name: str + type: type | None + default: Any + doc: str | None + + +def make_field_list(fields: list[Field]) -> nodes.Node: + node = addnodes.desc() + for f in fields: + f_node = addnodes.desc() + node.append(f_node) + signode = addnodes.desc_signature("", "") + f_node.append(signode) + signode += addnodes.desc_name(f.name, f.name) + if f.type is not None: + signode += addnodes.desc_annotation( + directives.unchanged, + "", + addnodes.desc_sig_punctuation("", ": "), + addnodes.desc_sig_space(), + nodes.Text(f.type.__name__ if isinstance(f.type, type) else str(f.type)), + ) + if f.default is not dataclasses.MISSING: + signode += addnodes.desc_annotation( + directives.unchanged, + "", + addnodes.desc_sig_punctuation("", " ="), + addnodes.desc_sig_space(), + nodes.Text(f.default), + ) + if f.doc is not None: + f_node.append(addnodes.desc_content("", nodes.Text(f.doc))) + + return node + + +def make_simple_field_list(fields: list[Field], *, defaults: bool = False) -> nodes.Node: + node = nodes.bullet_list() + for f in fields: + item = nodes.list_item() + node += item + p = nodes.paragraph() + p += nodes.literal(text=f.name) + text_pieces = [f.doc, f"(default: {f.default})" if f.default is not dataclasses.MISSING and defaults else None] + if f.doc is not None: + p += nodes.Text(": " + " ".join([piece for piece in text_pieces if piece])) + item += p + + return node + + def unindent(docstring: str) -> str: if not docstring: return "" @@ -139,12 +272,65 @@ def unindent(docstring: str) -> str: def line_indent(line: str) -> int | None: + """Determine the indent of a lines""" stripped = line.lstrip() if stripped: return len(line) - len(stripped) return None +def read_dataclass(c: type) -> list[Field]: + """Read the fields of a dataclass including docstrings for attributes.""" + docs = read_field_docstrings(c) + types = get_type_hints(c) + fields = dataclasses.fields(c) + return [Field(name=f.name, default=f.default, doc=docs.get(f.name), type=types.get(f.name)) for f in fields] + + +def read_field_docstrings(c: type) -> dict[str, str]: # noqa C903 + """Read field docstrings from a dataclass.""" + src = inspect.getsource(c) + indent = ((line_indent(src) or 0) + 4) * " " + + def find_line_start(src: str) -> str | None: + pos = src.find("\n" + indent) + return None if pos == -1 else src[pos + len(indent) + 1 :] + + def field_name(line: str) -> tuple[str, str | None]: + try: + name, rest = line.split(": ", 1) + except ValueError: + return line, None + return (rest, name) if name.isidentifier() else (line, None) + + def docstr(line: str) -> tuple[str, str | None]: + if not line.startswith('"""'): + return line, None + pos = line.find('"""', 3) + if pos == -1: + raise ValueError("Unterminated docstring found") + return line[pos + 3 :], line[3:pos] + + def tokenize(src: str) -> Iterable[tuple[str, str]]: + rest: str | None = src + f_name: str | None = None + while rest: + rest = find_line_start(rest) + if rest is None: + break + rest, new_f_name = field_name(rest) + if new_f_name is not None: + f_name = new_f_name + continue + if f_name is not None: + rest, docstring = docstr(rest) + if docstring is not None: + yield f_name, docstring + f_name = None + + return dict(tokenize(src)) + + @contextlib.contextmanager def mock_imports(whitelist: frozenset[str], extra_members: dict[str, list[str]] | None = None) -> Iterator[None]: """Add an import hook just after the builtin modules hook and the frozen module hook: diff --git a/docs/problems/airfoil.md b/docs/problems/airfoil.md index df423dd9..98b0b9d0 100644 --- a/docs/problems/airfoil.md +++ b/docs/problems/airfoil.md @@ -9,5 +9,165 @@ This issue does not occur on [Euler](https://docs.hpc.ethz.ch/#what-is-euler). On local machines, we recommend using [podman](https://podman.io/) or [docker](https://www.docker.com/) as the preferred runtime. ``` -``` {problem} airfoil +``` {problem:table} +:lead: Cashen Diniz @cashend ``` + +## Motivation + +The field of aerodynamics has always been a challenging testbed for engineering problems. In fact, many optimization and design methodologies were originally developed or perfected specifically for aerodynamic design applications [1]. Part of the reason for this is that even relatively simple aerodynamics problems can be complex, with slight changes in design parameters typically resulting in large changes in performance. In addition, the potential applications derived from solving these problems are quite practical, ranging from fixed-wing aircraft to hydrofoils used in ships, and wind turbine blades [2]. Here, we present Airfoil, a simple yet sufficiently realistic 2-dimensional aerodynamics benchmark problem. + +The airfoil problem presents a simple aerodynamic shape optimization routine based on Reynolds' averaged Navier-Stokes equations (RANS). In this problem, the solver attempts to indirectly morph the initial geometry to achieve a certain prescribed lift coefficient (which could correspond to a hypothetical loading requirement) while minimizing the amount of drag generated by the design. + +## Design space +The design space is represented as a tuple containing 192 2D points describing the airfoil coordinates and a scalar describing the rotation of the coordinates relative to the chord line needed to achieve a certain incoming direction of flow (the angle of attack, alpha): + +$$\text{Design Space} = \left\{\left(x, y\right)^{192}, \alpha \right\}$$ + +This specific coordinate parameterization (192) and rotational scalar design parameterization have been previously used in [3]. Another benchmark airfoil optimization problem, mentioned in [4] (Transonic RAE2822), was also used to internally validate the meshing and design parameterization. All training data, as well as the original and complete formulation of this problem can be found in [5]. + +## Objectives + +The objective is to minimize the coefficient of drag, $c_d$, and the optimization problem is defined as follows: + +$$\begin{aligned} +\min_{\Delta y_i,\alpha} \quad & c_d\\ +\text{s.t.} \quad & c_l = c_l^{\text{con}} \\ + & -0.025 \leq \Delta y_i \leq 0.025 \\ + & 0.0 \leq \alpha \leq 10.0 \\ + &\left( \frac{A}{A_{\text{init}}} \right)_{\text{min}} \leq \frac{A}{A_{\text{init}}} \leq 1.2 \\ +\end{aligned}$$ + +The terms used in the definition of the problem are described in Table 1. Note that some variables are defined relative to a required initial design input. + +**Table 1: Optimization Problem Parameters** + +| Category | Parameter | Quantity | Lower | Upper | Units | Description | +|----------|-----------|----------|-------|-------|-------|-------------| +| **Objective** | $c_d$ | 1 | - | - | Non-Dim./Counts | Coefficient of drag | +| **Variable** | $\Delta y_i$ | 20 | -0.025 | 0.025 | m | Change from initial FFD cage y value: $\Delta y_i = y_i - y_{\text{init}}$ | +| | $\alpha$ | 1 | 0.0 | 10.0 | Degrees | Angle of Attack | +| **Constraint** | $c_l = c_l^{\text{con}}$ | 1 | 0.0 | 0.0 | Non-Dim. | Coefficient of lift | +| | $\frac{A}{A_{\text{init}}}$ | 1 | $\left( \frac{A}{A_{\text{init}}} \right)_{\text{min}}$ | 1.20 | Non-Dim. | Area fraction; relative to initial | + +Instead of directly parameterizing all (192) coordinates of the airfoil, a smaller set of 20 control points is used. The collection of these control points forms what is known as a free-form deformation cage (FFD). Changes in the values of the FFD cage smoothly deform the underlying coordinates. In the airfoil problem, modifications to the geometry are parameterized by changes in the FFD control point y coordinates, $\Delta y_i$. See the EngiBench paper for more details. + +Note that, for simplicity, in this definition we have omitted certain constraints pertaining to thickness as well as those concerned with the shearing of the leading (front) and trailing (tail end) edges. + +## Conditions +The conditions are defined by the following parameters: + +```{problem:conditions} +``` + +Note that while it may be possible to constrain the design's area ahead of time, this is not typically possible for the prescribed coefficient of lift. In addition it may not be possible to achieve certain combinations of prescribed area and coefficient of lift constraints. + +## Simulator +All simulations used the open source and differentiable ADflow solver [6,7] as part of the MACH-Aero framework. ADflow was configured to run RANS simulations with the Spalart-Allmaras model for turbulence effects. Furthermore, within ADflow, the approximate Newton-Krylov method was used to improve convergence and robustness [8]. pyHyp, a hyperbolic mesh generator [9] was used to generate volume meshes automatically. For the optimization problem itself, we use the sequential least squares programming algorithm as implemented in the sparse optimization framework, pyOptSparse [10]. For geometry parameterization and deformation, module, we used the pyGeo and IDWarp frameworks [11,9]. + +You can install gcc and gfortran on your system with your package manager. +- On Ubuntu: `sudo apt-get install gcc gfortran` +- On MacOS: `brew install gcc gfortran` +- On Windows (WSL): `sudo apt install build-essential` + +## Dataset + +### v0 +The dataset linked to this problem is hosted on the [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/airfoil_v0). + +#### Fields +The dataset contains optimal design, conditions, objectives and these additional fields: +- `initial_design`: Design before the adjoint optimization. +- `cl_con_violation`: # Constraint violation for coefficient of lift. +- `area_ratio`: # Area ratio for given design. + +#### Creation Method +A dataset, originally described in [5], is integrated within our framework. The limits for the parameters in the data set are listed in Table 2. 1400 parameter combinations were chosen using Latin hypercube sampling (LHS) in the 4-dimensional parameter space. The training, testing and validation sets were randomly split into 748, 140, and 47 airfoil samples, respectively. Finally, Table 3 describes each variable in the data set available through the EngiBench API. + +**Table 2: Sampled Parameter Bounds** + +| Category | Parameter | Lower | Upper | Description | +|----------|-----------|-------|-------|-------------| +| **Flow Condition** | M | 0.4 | 0.9 | Mach number | +| | Re | 1E6 | 10E6 | Reynolds number | +| **Constraint** | $C_l^{\text{con}}$ | 0.5 | 1.2 | Coefficient of lift | +| | $\left( \frac{A}{A_{\text{init}}} \right)_{\text{min}}$ | 0.75 | 1.0 | Minimum area ratio | + +**Table 3: Airfoil dataset variables** + +| Category | Variable Name | Description | +|----------|---------------|-------------| +| **Flow Condition** | mach | Mach number | +| | reynolds | Reynolds number | +| **Constraint** | cl_target | Coefficient of lift (targeted constraint during optimization) | +| | area_ratio_min | Minimum area ratio | +| **Design Value** | area_initial | Area of the initial design | +| | cd | Coefficient of drag for the current design | +| | cl | Coefficient of lift for the selected design | +| | area_ratio | Area ratio of the selected design | + +## References +If you use this problem in your research, please cite the following paper: +```bibtex +@inproceedings{diniz2024optimizing, + author = {Diniz, C. and Fuge, M.}, + title = {Optimizing Diffusion to Diffuse Optimal Designs}, + booktitle = {AIAA SCITECH 2024 Forum}, + year = {2024}, + doi = {10.2514/6.2024-2013} +} +``` +[1] +Martins, J. R. R. A. and Ning, A. (2022). +Engineering Design Optimization. +Cambridge University Press, Cambridge, UK. + +[2] +Martins, J. (2022). +Aerodynamic design optimization: Challenges and perspectives. +Computers & Fluids, 239, 105391. + +[3] +Chen, W., Chiu, K., and Fuge, M. (2019). +Aerodynamic Design Optimization and Shape Exploration using Generative Adversarial Networks. +In AIAA Scitech 2019 Forum. + +[4] +He, X., Li, J., Mader, C. A., Yildirim, A., and Martins, J. R. R. A. (2019). +Robust aerodynamic shape optimization—from a circle to an airfoil. +Aerospace Science and Technology, 87, 48-61. + +[5] +Diniz, C. and Fuge, M. (2024). +Optimizing Diffusion to Diffuse Optimal Designs. +In AIAA SCITECH 2024 Forum. American Institute of Aeronautics and Astronautics. + +[6] +Mader, C. A., Kenway, G. K. W., Yildirim, A., and Martins, J. R. R. A. (2020). +ADflow—An open-source computational fluid dynamics solver for aerodynamic and multidisciplinary optimization. +Journal of Aerospace Information Systems. + +[7] +Kenway, G. K. W., Mader, C. A., He, P., and Martins, J. R. R. A. (2019). +Effective Adjoint Approaches for Computational Fluid Dynamics. +Progress in Aerospace Sciences, 110, 100542. + +[8] +Yildirim, A., Kenway, G. K. W., Mader, C. A., and Martins, J. R. R. A. (2019). +A Jacobian-free approximate Newton-Krylov startup strategy for RANS simulations. +Journal of Computational Physics, 397, 108741. + +[9] +Secco, N., Kenway, G. K. W., He, P., Mader, C. A., and Martins, J. R. R. A. (2021). +Efficient Mesh Generation and Deformation for Aerodynamic Shape Optimization. +AIAA Journal. + +[10] +Wu, N., Kenway, G., Mader, C. A., Jasa, J., and Martins, J. R. R. A. (2020). +pyOptSparse: A Python framework for large-scale constrained nonlinear optimization of sparse systems. +Journal of Open Source Software, 5(54), 2564. + +[11] +Hajdik, H. M., Yildirim, A., Wu, N., Brelje, B. J., Seraj, S., Mangano, M., Anibal, J. L., Jonsson, E., Adler, E. J., Mader, C. A., Kenway, G. K. W., and Martins, J. R. R. A. (2023). +pyGeo: A geometry package for multidisciplinary design optimization. +Journal of Open Source Software, 8(87), 5319. diff --git a/docs/problems/beams2d.md b/docs/problems/beams2d.md index b0acd185..e112ef28 100644 --- a/docs/problems/beams2d.md +++ b/docs/problems/beams2d.md @@ -1,4 +1,105 @@ # Beams2D -``` {problem} beams2d +``` {problem:table} +:lead: Arthur Drake @arthurdrake1 +``` + +## Motivation +The optimization of beam cross-sections is one of a fundamental problem in engineering, aiming to +maximize the structural stiffness under some applied force. This objective is usually formulated as +minimizing the compliance, which is the inverse of stiffness. In particular, TO frames the problem as +one of optimal material distribution, defining a grid of elements for which the material densities must +be determined on a scale from 0 to 1, where 1 represents the presence of material. After applying the +beam loads and other boundary conditions, designs are typically optimized using a gradient-based +approach with the help of the finite element method (FEM). While this is one of the simplest +TO applications, it is still a computationally expensive process requiring many iterations, opening the +door for faster approximation methods such as generative inverse design. +One of the most common beam types in TO is the Messerschmitt-Bölkow-Blohm (MBB) beam, +which is supported at the bottom-right and bottom-left corners, with a downward force applied +on the top-center. Given this symmetric configuration, one half of the design may be optimized +while representing the entire structure. We implement the MBB beam in ENGIBENCH for the most +accessible comparison to previous works in this domain. + +## Design Space +This problem simulates the right half-section of a MBB beam under bending. This half-beam is +subjected to a force at its top-left corner (corresponding to the top-center of the entire design) which +may also be shifted to the right to simulate different loading conditions. A roller support at the +bottom-right corner prevents vertical movement, and a symmetric boundary condition is enforced on +the left edge. The design space is an array of solid densities in `[0., 1.]` with a default size of +`(100, 50)` used by default, where `nelx = 100` and `nely = 50`. Internally, this is represented as +a flattened `(5000,)` array. Alternative shapes include `(50, 25)` for faster computation and `(200, 100)` +for higher-resolution results. Corresponding datasets for these three resolutions are provided. + +## Objectives +The goal is to optimize the distribution of solid material to minimize compliance +(equivalently, maximize stiffness) while satisfying constraints on material usage and minimum feature size. +Compliance is calculated as the sum of strain energy over the structure. + +The objectives are defined and indexed as follows: + +0. `c`: Compliance to minimize. + +## Conditions +The following input parameters define the problem conditions: + +```{problem:conditions} +``` + +## Simulator +Our simulation code is based on a Python adaptation of the popular 88-line topology optimization +code. It uses the more versatile density filtering approach in combination with a standard +Optimality Criteria (OC) optimization method. Two primary sensitivity matrices, one with respect +to compliance (`dc`) and the other with respect to volume fraction (`dv`), are continuously updated +and used to calculate a given design's compliance value. We have also ensured that during the +required Lagrange multiplier search within OC, the inner optimization loop terminates if the absolute +difference upper and lower bounds diminishes to a value smaller than machine precision. This +prevents the code from becoming stuck at this point, which we observed in some warm-starting +instances with noisy initial designs. + +Compliance `c` is calculated using: +```python +c = ((Emin + xPrint**penal * (Emax - Emin)) * ce).sum() +``` + +where `xPrint` is the current true density field, `penal` is the penalization factor (e.g., 3.0), +and `ce` is the element-wise strain energy density. + +## Dataset +This problem offers multiple datasets for various sizes of `nelx` and `nely`. Each dataset includes +columns for the optimal design, all conditions listed above, and the corresponding objective values. +For advanced usage, we also provide a column containing the optimization history. The datasets have +been generated by sampling conditions over a structured grid for various problem sizes. +Three datasets are available on the [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab). +They correspond to resolutions of $50 \times 25$, $100 \times 50$ (default), and $200 \times 100$. + +### v0 + +#### Fields +Each dataset contains: +- Optimized beam structures, +- The corresponding condition parameters, +- Objective values (compliance), +- Full optimization histories (for advanced use). + +#### Creation Method +Datasets were generated by uniformly sampling the condition space. The resolutions used are: +- `(50, 25)` +- `(100, 50)` +- `(200, 100)` + +A more comprehensive description of the creation method can be found in the [README](https://github.com/IDEALLab/EngiBench/tree/main/engibench/problems/beams2d). + +## Citation +If you use this problem in your research, please cite the following paper: +``` +@article{andreassen2011efficient, + title={Efficient topology optimization in MATLAB using 88 lines of code}, + author={Andreassen, Erik and Clausen, Anders and Schevenels, Mattias and Lazarov, Boyan S and Sigmund, Ole}, + journal={Structural and Multidisciplinary Optimization}, + volume={43}, + number={1}, + pages={1--16}, + year={2011}, + publisher={Springer} +} ``` diff --git a/docs/problems/heatconduction2d.md b/docs/problems/heatconduction2d.md index 0f288751..d6b99b8a 100644 --- a/docs/problems/heatconduction2d.md +++ b/docs/problems/heatconduction2d.md @@ -1,4 +1,40 @@ # Heat Conduction 2D -``` {problem} heatconduction2d +``` {problem:table} +:lead: Milad Habibi @MIladHB ``` + +## Design space +These problems are a specific subset of topology optimization problems aimed at minimizing thermal compliance within a unit square (2D) subject to: a constraint on the volume of highly conductive material used, and given boundary conditions, +particularly the location of the adiabatic region. The adiabatic region refers to a symmetric length on the bottom side of the 2D problem space . The design space for the 2D problem consists of a 2D array representing solid densities, +which is parametrized by resolution, that is, 'DesignSpace = [0,1] By default, a $101 \times 101$ space is used for the 2D problem. + +## Objectives +The objective is defined and indexed as follows: + +0. `c`: Thermal compliance coefficient to minimize. + +## Conditions +The conditions are defined by the following parameters: + +```{problem:conditions} +``` + +## Simulator +The simulator is a docker container with the dolfin-adjoint software that computes the thermal compliance of the design. +We convert use intermediary files to convert from and to the simulator that is run from a Docker image. + +## Dataset +The dataset has been generated the dolfin-adjoint software. It is hosted on the [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/heat_conduction_2d_v0). + +### v0 + +#### Fields +The dataset only contains conditions and optimal designs (no objective). + +#### Creation Method +The creation method for the dataset is specified in the reference paper. + +## References +If you use this problem in your research, please cite the following paper: +Milad Habibi, Jun Wang, and Mark Fuge, "When Is it Actually Worth Learning Inverse Design?" in IDETC 2023. doi: https://doi.org/10.1115/DETC2023-116678 diff --git a/docs/problems/heatconduction3d.md b/docs/problems/heatconduction3d.md index f05844b5..91ab0207 100644 --- a/docs/problems/heatconduction3d.md +++ b/docs/problems/heatconduction3d.md @@ -1,4 +1,46 @@ # Heat Conduction 3D -``` {problem} heatconduction3d +``` {problem:table} +:lead: Milad Habibi @MIladHB ``` + +## Motivation +Heat conduction problems serve as fundamental benchmarks for the development and evaluation of design optimization methods, with applications ranging from thermal management in electronic devices to insulation systems and +heat exchangers in industrial applications. As thermal management has become critical in fields such as aerospace, automotive, and consumer electronics, both industry and academia have shown growing interest in advanced thermal +design systems. In response to this demand, topology optimization has become popular as a powerful approach for improving heat dissipation while minimizing material usage. In addition, the development of additive manufacturing +technologies has made the complex geometries produced by topology optimization more feasible to fabricate in real-world applications. + +## Design space +These problems are a specific subset of topology optimization problems aimed at minimizing thermal compliance within a unit cube (3D), subject to: a constraint on the volume of highly conductive material used, and given boundary conditions, +particularly the location of the adiabatic region. The adiabatic region refers to a prescribed symmetric area on the bottom surface of the 3D problem space. The 3D design space is represented as a 3D tensor of densities. By default, +a $51 \times 51 \times 51$ space is used for the 3D problem. + +## Objectives +The objective is defined and indexed as follows: + +0. `c`: Thermal compliance coefficient to minimize. + +## Conditions +The conditions are defined by the following parameters: + +```{problem:conditions} +``` + +## Simulator +The simulator is a docker container with the dolfin-adjoint software that computes the thermal compliance of the design. +We convert use intermediary files to convert from and to the simulator that is run from a Docker image. + +## Dataset +The dataset has been generated the dolfin-adjoint software. It is hosted on the [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/heat_conduction_3d_v0). + +### v0 + +#### Fields +The dataset only contains conditions and optimal designs (no objective). + +#### Creation Method +The creation method for the dataset is specified in the reference paper. + +## References +If you use this problem in your research, please cite the following paper: +Habibi, Milad, Shai Bernard, Jun Wang, and Mark Fuge, "Mean squared error may lead you astray when optimizing your inverse design methods" in JMD 2025. doi: https://doi.org/10.1115/1.4066102 diff --git a/docs/problems/photonics2d.md b/docs/problems/photonics2d.md index db1637b3..93525f0a 100644 --- a/docs/problems/photonics2d.md +++ b/docs/problems/photonics2d.md @@ -1,4 +1,150 @@ # Photonics2D -``` {problem} photonics2d +``` {problem:table} +:lead: Mark Fuge @markfuge +``` + +## Motivation +The optimization of photonic circuits in general, and multiplexers in particular, was +one of the initial and most widely studied problems in the inverse design of +electromagnetic/optical devices. In part, this is because multiplexer +devices have several interesting properties that make them more difficult to create generative +models of, compared to other problems in EngiBench. This includes the fact that, due to the wave +properties of the physical phenomena there are usually multiple solutions with equivalent or +similar performance, which results from shifting or inverting the phase profile of the +electromagnetic wave. This adds complexity to the generative model in that the solution may not +have a single unique global minimum. Another motivating factor for including this problem is the +complexity of the structures/designs themselves: unlike in structural or thermal compliance +problems, which lead to connected structures, the photonics solutions often involve several +disconnected elements whose relative position and spacing is governed by the specific wavelengths +it needs to demultiplex. This is a difficult prediction and generation task, compared to, e.g., +generating a connected beam structure. Thus, it acts as a good counterpoint to add to the library +and provides a mechanism to benchmark generative algorithms that can perform well on both +connected and disconnected design topologies. + +## Design space +This problem simulates a wavelength demultiplexer where the optimized device will direct an +electromagnetic wave to one of two possible output ports depending on the wavelength/frequency of +the incoming wave. Specifically, the demultiplexer targets two specific wavelengths (referred to +$\lambda_1$ and $\lambda_2$ in the library), and the performance of the device is how well it can +bend or direct the energy toward two specific locations in the device, as measured by how much of +the electric field of each wavelength overlaps with the desired output port locations. The design +space consists of a 2D array representing the presence of either a high or low permittivity, +parameterized by `nelx` and `nely`, i.e., $\text{design_space} = [0,1]^{\text{nelx}\times \text{nely}}$. +By default, the library uses a $120 \times 120$ space, however, this can be modified to non-square +design spaces by the user. Specifically, we use a 2D tensor `rho` (num_elems_x, num_elems_y) with +values in [0, 1], representing material density. This is stored as `design_space` (gymnasium.spaces.Box). + +## Objectives +The main objective is to maximize the overlap of the electric field of the simulated wavelength at +the target output location, with an optional penalty for the amount of material used (this penalty +weight is set to a small default value ($1e^{-2}$) for consistency, but can be altered for advanced +usage): +0. `total_overlap`: Objective to maximize, defined as + `overlap1 * overlap2`. Higher is better. This corresponds to + the overlap in the target electrical fields with the desired demultiplexing locations. + Note that bot `simulate` and `optimize` subtract a small material penalty + (`total_overlap - penalty`) to avoid multiple equivalent local optima, but this penalty + is small relative to the overlap objective. + +## Conditions +These are designed as user-configurable parameters that alter the problem definition. The +conditions consist of the two input wavelengths to be demultiplexed -- +$\lambda_1$ and $\lambda_2$, as well as a desired `blur_radius` ($r_{blur}$) parameter, +which blurs (using a circular convolution) the pixelized design field for a chosen number of +integer pixels\textemdash this blurring essentially creates a penalty on the minimum feature size +of the design. The size of the device -- expressed as `nelx` and `nely` -- is also adjustable, +and could be viewed as a possible condition for multi-resolution problems, but in practice, as with +`Beams2D`, this is built into the problem definition since it produces a different dataset. + +Default problem parameters that can be overridden via the `config` dict: + +```{problem:conditions} +:defaults: true +``` + +In practice, for the dataset loading, we will keep `num_elems_x` and `num_elems_y`to set +values for each dataset, such that different resolutions correspond to different +independent datasets. + +## Optimization Parameters +Note: These are advanced parameters that alter the optimization process -- +we do not recommend changing these if you are only using the library for benchmarking, +as it could make results less reproducible across papers using this problem.) +- `num_optimization_steps`: Total number of optimization steps (default: 300). +- `step_size`: Adam optimizer step size (default: 1e-1). +- `penalty_weight`: Weight for the L2 penalty term (default: 1e-2). Larger values reduce + unnecessary material, but may lead to worse performance if too large. +- `eta`: Projection center parameter (default: 0.5). There is little reason to change this. +- `N_proj`: Number of projection applications (default: 1). Increasing this can help make + the design more binary. +- `N_blur`: Number of blur applications (default: 1). Increasing this smooths the design more. +- `initial_beta`: Initial beta for the optimization continuation scheme (default: 1.0). +- `save_frame_interval`: Interval for saving intermediate design frames during optimization. + If > 0, saves a frame every `save_frame_interval` iterations + to the `opt_frames/` directory. Default is 0 (disabled). + +## Internal Constants +Note: These are not typically changed by users, but provided here for technical reference +- `dl`: Spatial resolution (meters) (default: 40e-9). +- `Npml`: Number of PML cells (default: 20). +- `epsr_min`: Minimum relative permittivity (default: 1.0). +- `epsr_max`: Maximum relative permittivity (default: 12.0). +- `space_slice`: Extra space for source/probe slices (pixels) (default: 8). + +## Simulator +The simulation uses the `ceviche` library's Finite Difference Frequency Domain (FDFD) +solver (`fdfd_ez`). Optimization uses `ceviche.optimizers.adam_optimize` with +gradients computed via automatic differentiation (`autograd`). +The simulation code uses the \texttt{ceviche} library and specifically, +the wave demultiplexer demonstration case provided by the +[library authors](https://github.com/fancompute/workshop-invdesign/blob/master/04_Invdes_wdm_scheduling.ipynb} +based on their related publication, which uses a similar formalism to an earlier demultiplexer +paper by Piggott (2015). The optimization method is first-order and uses the Adam optimizer. Beyond +the baseline implementation already available via `ceviche`, we implemented a polynomial $\beta$ continuation +scheme that performed more reliably than the step-wise continuation scheme used in the original implementation, +and EngiBench also possesses the ability to change the starting and ending continuation values, for future +research cases where one wishes to estimate or optimize the continuation schedule themselves. Other than these +changes, the implementation of this problem is as consistent as possible with that of the original `ceviche` library. + +## Dataset +This problem offers a single datasets of `nelx`=120 and `nely`=120, although various sizes of `nelx` and `nely` +could be generated from the library if desired. The dataset includes columns for the optimal design, all conditions +listed above, and the corresponding objective value history as the optimizer optimized toward the optimal design +provided in the dataset. The dataset was generated by sampling by sampling at random $\lambda_1$, $\lambda_2$ and +$r_{blur}$ over the conditions mentioned above. The dataset is available on the +[Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/photonics_2d_120_120_v0). + +### v0 + +#### Fields +Each dataset contains: +- `lambda1`: The first input wavelength in μm. +- `lambda2`: The second input wavelength in μm. +- `blur_radius`: Radius for the density blurring filter (pixels). +- `optimal_design`: The optimal design density array (shape num_elems_x, num_elems_y). +- `optimization_history`: A list of objective values from the optimization process (field overlap minus penalty, where higher is better) -- This is for advanced use. + +#### Creation Method +To generate a dataset for training, we generate (randomly, uniformly) swept over the following parameters: +- $\lambda_1 \in [0.5\mu m, 1.25\mu m]$ = `lambda1` = `rng.uniform(low=0.5, high=1.25, size=20)` -- This corresponds roughly to a portion of the visible spectrum up to near-infrared. +- $\lambda_2 \in [0.75\mu m, 1.5\mu m]$ = `lambda2` = `rng.uniform(low=0.75, high=1.5, size=20)` -- This corresponds roughly to a portion of the visible spectrum up to near-infrared. +- $r_{blur}$ = `blur_radius` = `range(0, 5)` + +## Citation +This problem is directly refactored from the Ceviche Library: +https://github.com/fancompute/ceviche +and if you use this problem your experiments, you can use the citation below +provided by the original library authors: +``` +@article{hughes2019forward, + title={Forward-Mode Differentiation of Maxwell's Equations}, + author={Hughes, Tyler W and Williamson, Ian AD and Minkov, Momchil and Fan, Shanhui}, + journal={ACS Photonics}, + volume={6}, + number={11}, + pages={3010--3016}, + year={2019}, + publisher={ACS Publications} +} ``` diff --git a/docs/problems/power_electronics.md b/docs/problems/power_electronics.md index c324769a..b7c6c65a 100644 --- a/docs/problems/power_electronics.md +++ b/docs/problems/power_electronics.md @@ -1,4 +1,81 @@ # PowerElectronics -``` {problem} power_electronics +``` {problem:table} +:lead: Xuliang Dong @ liangXD523 ``` + +## Motivation +Optimizing circuit parameters is a critical aspect of circuit design but remains challenging, particularly for power converter circuits that contain diodes and switches, +which introduce significant nonlinearity and discontinuity. These characteristics make key objectives such as *DcGain* and *Voltage Ripple* highly sensitive to even small parameter variations. + +Because the circuit simulator NgSpice operates as a black box and is non-differentiable, gradient-based optimization methods are not suitable. +Bayesian optimization is commonly employed for parameter tuning, while surrogate models offer a promising alternative. +Even under the constraint of a fixed topology, optimizing circuit parameters to minimize the objectives remains a difficult problem for surrogate models. + +NgSpice applies transient analysis by formulating the system as a set of differential equations based on Kirchhoff's laws. +These equations are discretized using numerical integration methods such as the Backward Euler or trapezoidal rule and solved iteratively at each time step to compute performance metrics. +To ensure stable simulations, a specific on-off switching pattern is chosen for the circuit. +Despite this simplification, determining the optimal parameter values remains highly challenging. + +## Design Space +The design space for this problem is represented as a 10-dimensional bounded box, where each dimension corresponds to a specific circuit parameter. These parameters include values for capacitors, inductors, and a shared duty cycle for all switches. Each design can be expressed as a vector **x** of the form: + +$$ +x = \begin{bmatrix} C_1,\dots,C_6,L_1,L_2,L_3,T_1 \end{bmatrix}^{\top} \in \mathcal{X}, +\quad +\mathcal{X} = [1\text{e}{-6}, 2\text{e}{-5}]^6 \times [1\text{e}{-6}, 1\text{e}{-3}]^3 \times [0.1, 0.9] +$$ + +Here, $C_1,\dots,C_6$ are the capacitance values (in Farads), $L_1,L_2,L_3$ are the inductance values (in Henries), and $T_1$ is the duty cycle shared across all 5 switches. The duty cycle $T_1$ denotes the fraction of time during which the switches are in the “on” state and governs a periodic on-off pattern repeated at high frequency throughout the simulation. + +## Objectives +The simulation outputs two scalar values: *DcGain* and *Voltage Ripple*. The former represents the ratio of load to input voltage and should ideally approximate a predefined constant, such as $0.25$, as closely as possible. Meanwhile, the latter quantifies the voltage fluctuation at the load. + +**DcGain objective:** + +$$ +\min_{\mathbf{x} \in \mathcal{X}} \; \bigg|\frac{\overline{V_{load}(t)}}{V_{source}} - 0.25\bigg| += \bigg|\frac{1}{V_{source}} \cdot \frac{1}{T} \sum_{i=1}^{N-1} \frac{V_{load}(t_{i+1}) + V_{load}(t_i)}{2} \cdot (t_{i+1} - t_i) - 0.25\bigg| +$$ + +where $\overline{V_{load}(t)}$ is the average load voltage, $V_{source} = 1000$ volts, and $T = t_N - t_1$ is the simulation duration. + +**Voltage Ripple objective:** + +$$ +\min_{\mathbf{x} \in \mathcal{X}} \; \text{Voltage Ripple} += \frac{V_{pp}(t)}{\overline{V_{load}(t)}} += \frac{\max_{i \in [1, N]} V_{load}(t_i) - \min_{i \in [1, N]} V_{load}(t_i)}{\overline{V_{load}(t)}} +$$ + +where $V_{pp}$ is the peak-to-peak load voltage calculated during transient analysis. + +## Conditions +This problem does not include environmental or operational conditions as part of its input specification. Unlike other domains where the simulation setup may vary based on conditions (e.g., load configurations or external temperatures), the circuit is simulated under fixed source voltage and switching behavior. As a result, the design optimization task focuses solely on tuning internal circuit parameters, with no external conditions to vary. More complex variants of this problem — involving multiple topologies or variable source voltages — may be considered in future releases. + +## Simulator +The simulator is ngSpice circuit simulator. You can download it based on your operating system: +- Windows: [https://sourceforge.net/projects/ngspice/files/ng-spice-rework/45.2/](https://sourceforge.net/projects/ngspice/files/ng-spice-rework/45.2/) +- MacOS: `brew install ngspice` +- Linux: `sudo apt-get install ngspice` + +## Dataset +The dataset linked to this problem is hosted on the [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/power_electronics). + +### v0 + +#### Fields +The dataset contains 3 fields: +- `initial_design`: The 20-dimensional design variable defined above. +- `DcGain`: The ratio of load vs. input voltage. +- `Voltage_Ripple`: The fluctuation of voltage on the load `R0`. + +#### Creation Method +We created this dataset in 3 parts. All the 3 parts are simulated with {`GS0_L1`, `GS1_L1`, `GS2_L1`, `GS3_L1`, `GS4_L1`} = {1, 0, 0, 1, 1} and {`GS0_L2`, `GS1_L2`, `GS2_L2`, `GS3_L2`, `GS4_L2`} = {1, 0, 1, 1, 0}. +Here are the 3 parts: +1. 6 capacitors and 3 inductors only take their min and max values. `T1` ranges {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}. There are 2^6 * 2^3 * 9 = 4608 samples. +2. Random sample 4608 points in the 6 + 3 + 1 = 10 dimensional space. Min and max values in each dimension will not be sampled. +3. Latin hypercube sample 4608 points in the 6 + 3 + 1 = 10 dimensional space. Each dimension is split into 10 intervals. Min and max values in each dimension will not be sampled. + +## References +If you use this problem in your research, please cite the following paper: diff --git a/docs/problems/thermoelastic2d.md b/docs/problems/thermoelastic2d.md index 0a30b433..ebdc8299 100644 --- a/docs/problems/thermoelastic2d.md +++ b/docs/problems/thermoelastic2d.md @@ -1,4 +1,54 @@ # Thermoelastic2D -``` {problem} thermoelastic2d +```{problem:table} +:lead: Gabriel Apaza @gapaza ``` + +## Motivation +As articulated in their respective sections, both the Beams2D and HeatConduction2D problems found in the EngiBench library are fundamental engineering design problems that have historically served as benchmarks for the development and testing of optimization methods. +While their relevance is supported by needs in real engineering design scenarios (aerospace, automotive, consumer electronics, etc...), their mono-domain nature ignores the reality that coupling between domains exists, and should be accounted for in scenarios where performance in one domain significantly impacts performance in another. +To address this distinction, a multi-physics topology optimization problem is developed that captures the coupling between structural and thermal domains. + +## Design space +This multi-physics topology optimization problem is governed by linear elasticity and steady-state heat conduction with a one-way coupling from the thermal domain to the elastic domain. +The problem is defined over a square 2D domain, where load elements and support elements are placed along the boundary to define a unique elastic condition. +Similarly, heatsink elements are placed along the boundary to define a unique thermal condition. +The design space is then defined by a 2D array representing density values (parameterized by DesignSpace = [0,1]^{`nelx` x `nely`}, where nelx and `nely` denote the x and y dimensions). + +## Objectives +The objective of this problem is to minimize total compliance C under a volume fraction constraint V by placing a thermally conductive material. +Total compliance is defined as the sum of thermal compliance and structural compliance. + +## Conditions + +```{problem:conditions} +``` + +## Simulator +The simulation code is based on a Python adaptation of the popular 88-line topology optimization code, modified to handle the thermal domain in addition to thermal-elastic coupling. +Optimization is conducted by reformulating the integer optimization problem as a continuous one (leveraging a SIMP approach), where a density filtering approach is used to prevent checkerboard-like artifacts. +The optimization process itself operates by calculating the sensitivities of the design variables with respect to total compliance (done efficiently using the Adjoint method), calculating the sensitivities of the design variables with respect to the constraint value, and then updating the design variables by solving a convex-linear subproblem and taking a small step (using the method of moving asymptotes). +The optimization loop terminates when either an upper bound of the number of iterations has been reached or if the magnitude of the gradient update is below some threshold. + +## Dataset +The dataset linked to this problem is on huggingface [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/thermoelastic_2d_v0). +This dataset contains a set of 1000 optimized thermoelastic designs in a 64x64 domain, where each design is optimized for a unique set of conditions. +Each datapoint's conditions are randomly generated by arbitrarily placing: a single loaded element along the bottom boundary, two fixed elements (fixed in both the x and y direction) along the left and top boundary, and heatsink elements along the right boundary. +Furthermore, values for the volume fraction constraint are randomly selected in the range $[0.2, 0.5]$. + +Relevant datapoint fields include: +- `optimal_design`: An optimized design for the set of boundary conditions +- `strain`: The strain field for the initial uniform design +- `vm_stress`: The von Mises stress field for the initial uniform design +- `fixed_elements`: Encodes a binary NxN matrix of the structurally fixed elements in the domain. +- `force_elements_x`: Encodes a binary NxN matrix specifying elements that have a structural load in the x-direction. +- `force_elements_y`: Encodes a binary NxN matrix specifying elements that have a structural load in the y-direction. +- `heatsink_elements`: Encodes a binary NxN matrix specifying elements that have a heat sink. +- `volume_fraction`: The volume fraction value of the optimized design +- `structural_compliance`: The structural compliance of the optimized design +- `thermal_compliance`: The thermal compliance of the optimized design +- `nelx`: The number of elements in the x-direction +- `nely`: The number of elements in the y-direction +- `volfrac`: The volume fraction target of the optimized design +- `rmin`: The filter size used in the optimization routine +- `weight`: The domain weighting used in the optimization routine diff --git a/docs/tutorials/new_problem.md b/docs/tutorials/new_problem.md index d048ec00..a909454b 100644 --- a/docs/tutorials/new_problem.md +++ b/docs/tutorials/new_problem.md @@ -54,15 +54,44 @@ In general, follow the `beams2d/` example. 1. Install necessary documentation tools: `pip install ".[doc]"`. 2. If it is a new problem family, add a new `.md` file in [docs/problems/](source:docs/problems/) following the existing structure and add your problem family in the `toctree` of [docs/problems/index.md](source:docs/problems/index.md). -3. Add a problem markdown file to the `toctree` in `docs/problems/new_problem.md`. In the md file, use EngiBench's own `problem` directive: +3. Add a problem markdown file to the `toctree` in `docs/problems/new_problem.md`. In the md file, use EngiBench's own `problem:table` and `problem:conditions` directives in the docs of your new problem: ``````md # Your Problem - ``` {problem} new_problem + ``` {problem:table} + :lead: Chuck Norris @chucknorris ``` + + ... + + ## Conditions + + ``` {problem:conditions} + ``` + + ... `````` - Here, `new_problem` must match the name of the top level module where your problem class is defined. + **`problem:table`**: This directive extracts metadata from a problem + and inserts a table filled with the metadata. + By default, the directive will try to import the problem + `engibench.problems.`, where `` is the filename + (without `.md` extension) of the markdown file the directive is used. + + Options (optional): + * `:problem_id:` override ``, + * `:lead:` Add a row "Lead" to the table, containing the specified value. + If the value ends with `@username`, a link to `https://github.com/username` will be inserted. + + **`problem:conditions`**: This directive lists the conditions + extracted from a problem as in the "Conditions" row produced by the + `problem:table` directive. The `` is determined the same + way as in `problem:table`. + + Options: + * `:problem_id:` override ``, + * `:defaults:` include default values in the list of conditions + Here, `new_problem/__init__.py` is crucial as it makes the problem class discoverable to the `problem` directive by the reexport `from engibench.problems.new_problem.v0 import NewProblem`. 4. Add an image (result of `problem.render(design)`) in `docs/_static/img/problems`. The file's name should be `.png`, with your problem module as in the point above. diff --git a/engibench/problems/airfoil/v0.py b/engibench/problems/airfoil/v0.py index 6cd3916c..303270bd 100644 --- a/engibench/problems/airfoil/v0.py +++ b/engibench/problems/airfoil/v0.py @@ -83,172 +83,7 @@ def does_not_self_intersect(design: DesignType) -> None: class Airfoil(Problem[DesignType]): r"""Airfoil 2D shape optimization problem. - ## Problem Description This problem simulates the performance of an airfoil in a 2D environment. An airfoil is represented by a set of 192 points that define its shape. The performance is evaluated by the [MACH-Aero](https://mdolab-mach-aero.readthedocs-hosted.com/en/latest/) simulator that computes the lift and drag coefficients of the airfoil. - - ## Motivation - - The field of aerodynamics has always been a challenging testbed for engineering problems. In fact, many optimization and design methodologies were originally developed or perfected specifically for aerodynamic design applications [1]. Part of the reason for this is that even relatively simple aerodynamics problems can be complex, with slight changes in design parameters typically resulting in large changes in performance. In addition, the potential applications derived from solving these problems are quite practical, ranging from fixed-wing aircraft to hydrofoils used in ships, and wind turbine blades [2]. Here, we present Airfoil, a simple yet sufficiently realistic 2-dimensional aerodynamics benchmark problem. - - The airfoil problem presents a simple aerodynamic shape optimization routine based on Reynolds' averaged Navier-Stokes equations (RANS). In this problem, the solver attempts to indirectly morph the initial geometry to achieve a certain prescribed lift coefficient (which could correspond to a hypothetical loading requirement) while minimizing the amount of drag generated by the design. - - ## Design space - The design space is represented as a tuple containing 192 2D points describing the airfoil coordinates and a scalar describing the rotation of the coordinates relative to the chord line needed to achieve a certain incoming direction of flow (the angle of attack, alpha): - - $$\text{Design Space} = \left\{\left(x, y\right)^{192}, \alpha \right\}$$ - - This specific coordinate parameterization (192) and rotational scalar design parameterization have been previously used in [3]. Another benchmark airfoil optimization problem, mentioned in [4] (Transonic RAE2822), was also used to internally validate the meshing and design parameterization. All training data, as well as the original and complete formulation of this problem can be found in [5]. - - ## Objectives - - The objective is to minimize the coefficient of drag, $c_d$, and the optimization problem is defined as follows: - - $$\begin{aligned} - \min_{\Delta y_i,\alpha} \quad & c_d\\ - \text{s.t.} \quad & c_l = c_l^{\text{con}} \\ - & -0.025 \leq \Delta y_i \leq 0.025 \\ - & 0.0 \leq \alpha \leq 10.0 \\ - &\left( \frac{A}{A_{\text{init}}} \right)_{\text{min}} \leq \frac{A}{A_{\text{init}}} \leq 1.2 \\ - \end{aligned}$$ - - The terms used in the definition of the problem are described in Table 1. Note that some variables are defined relative to a required initial design input. - - **Table 1: Optimization Problem Parameters** - - | Category | Parameter | Quantity | Lower | Upper | Units | Description | - |----------|-----------|----------|-------|-------|-------|-------------| - | **Objective** | $c_d$ | 1 | - | - | Non-Dim./Counts | Coefficient of drag | - | **Variable** | $\Delta y_i$ | 20 | -0.025 | 0.025 | m | Change from initial FFD cage y value: $\Delta y_i = y_i - y_{\text{init}}$ | - | | $\alpha$ | 1 | 0.0 | 10.0 | Degrees | Angle of Attack | - | **Constraint** | $c_l = c_l^{\text{con}}$ | 1 | 0.0 | 0.0 | Non-Dim. | Coefficient of lift | - | | $\frac{A}{A_{\text{init}}}$ | 1 | $\left( \frac{A}{A_{\text{init}}} \right)_{\text{min}}$ | 1.20 | Non-Dim. | Area fraction; relative to initial | - - Instead of directly parameterizing all (192) coordinates of the airfoil, a smaller set of 20 control points is used. The collection of these control points forms what is known as a free-form deformation cage (FFD). Changes in the values of the FFD cage smoothly deform the underlying coordinates. In the airfoil problem, modifications to the geometry are parameterized by changes in the FFD control point y coordinates, $\Delta y_i$. See the EngiBench paper for more details. - - Note that, for simplicity, in this definition we have omitted certain constraints pertaining to thickness as well as those concerned with the shearing of the leading (front) and trailing (tail end) edges. - - ## Conditions - The conditions are defined by the following parameters: - - `mach`: Mach number. - - `reynolds`: Reynolds number. - - `area_ratio_min`: Minimum area ratio (ratio relative to initial area) constraint. - - `area_initial`: Initial area. - - `cl_target`: Target lift coefficient to satisfy equality constraint. - - Note that while it may be possible to constrain the design's area ahead of time, this is not typically possible for the prescribed coefficient of lift. In addition it may not be possible to achieve certain combinations of prescribed area and coefficient of lift constraints. - - ## Simulator - All simulations used the open source and differentiable ADflow solver [6,7] as part of the MACH-Aero framework. ADflow was configured to run RANS simulations with the Spalart-Allmaras model for turbulence effects. Furthermore, within ADflow, the approximate Newton-Krylov method was used to improve convergence and robustness [8]. pyHyp, a hyperbolic mesh generator [9] was used to generate volume meshes automatically. For the optimization problem itself, we use the sequential least squares programming algorithm as implemented in the sparse optimization framework, pyOptSparse [10]. For geometry parameterization and deformation, module, we used the pyGeo and IDWarp frameworks [11,9]. - - You can install gcc and gfortran on your system with your package manager. - - On Ubuntu: `sudo apt-get install gcc gfortran` - - On MacOS: `brew install gcc gfortran` - - On Windows (WSL): `sudo apt install build-essential` - - ## Dataset - - ### v0 - The dataset linked to this problem is hosted on the [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/airfoil_v0). - - #### Fields - The dataset contains optimal design, conditions, objectives and these additional fields: - - `initial_design`: Design before the adjoint optimization. - - `cl_con_violation`: # Constraint violation for coefficient of lift. - - `area_ratio`: # Area ratio for given design. - - #### Creation Method - A dataset, originally described in [5], is integrated within our framework. The limits for the parameters in the data set are listed in Table 2. 1400 parameter combinations were chosen using Latin hypercube sampling (LHS) in the 4-dimensional parameter space. The training, testing and validation sets were randomly split into 748, 140, and 47 airfoil samples, respectively. Finally, Table 3 describes each variable in the data set available through the EngiBench API. - - **Table 2: Sampled Parameter Bounds** - - | Category | Parameter | Lower | Upper | Description | - |----------|-----------|-------|-------|-------------| - | **Flow Condition** | M | 0.4 | 0.9 | Mach number | - | | Re | 1E6 | 10E6 | Reynolds number | - | **Constraint** | $C_l^{\text{con}}$ | 0.5 | 1.2 | Coefficient of lift | - | | $\left( \frac{A}{A_{\text{init}}} \right)_{\text{min}}$ | 0.75 | 1.0 | Minimum area ratio | - - **Table 3: Airfoil dataset variables** - - | Category | Variable Name | Description | - |----------|---------------|-------------| - | **Flow Condition** | mach | Mach number | - | | reynolds | Reynolds number | - | **Constraint** | cl_target | Coefficient of lift (targeted constraint during optimization) | - | | area_ratio_min | Minimum area ratio | - | **Design Value** | area_initial | Area of the initial design | - | | cd | Coefficient of drag for the current design | - | | cl | Coefficient of lift for the selected design | - | | area_ratio | Area ratio of the selected design | - - ## References - If you use this problem in your research, please cite the following paper: - ```bibtex - @inproceedings{diniz2024optimizing, - author = {Diniz, C. and Fuge, M.}, - title = {Optimizing Diffusion to Diffuse Optimal Designs}, - booktitle = {AIAA SCITECH 2024 Forum}, - year = {2024}, - doi = {10.2514/6.2024-2013} - } - ``` - [1] - Martins, J. R. R. A. and Ning, A. (2022). - Engineering Design Optimization. - Cambridge University Press, Cambridge, UK. - - [2] - Martins, J. (2022). - Aerodynamic design optimization: Challenges and perspectives. - Computers & Fluids, 239, 105391. - - [3] - Chen, W., Chiu, K., and Fuge, M. (2019). - Aerodynamic Design Optimization and Shape Exploration using Generative Adversarial Networks. - In AIAA Scitech 2019 Forum. - - [4] - He, X., Li, J., Mader, C. A., Yildirim, A., and Martins, J. R. R. A. (2019). - Robust aerodynamic shape optimization—from a circle to an airfoil. - Aerospace Science and Technology, 87, 48-61. - - [5] - Diniz, C. and Fuge, M. (2024). - Optimizing Diffusion to Diffuse Optimal Designs. - In AIAA SCITECH 2024 Forum. American Institute of Aeronautics and Astronautics. - - [6] - Mader, C. A., Kenway, G. K. W., Yildirim, A., and Martins, J. R. R. A. (2020). - ADflow—An open-source computational fluid dynamics solver for aerodynamic and multidisciplinary optimization. - Journal of Aerospace Information Systems. - - [7] - Kenway, G. K. W., Mader, C. A., He, P., and Martins, J. R. R. A. (2019). - Effective Adjoint Approaches for Computational Fluid Dynamics. - Progress in Aerospace Sciences, 110, 100542. - - [8] - Yildirim, A., Kenway, G. K. W., Mader, C. A., and Martins, J. R. R. A. (2019). - A Jacobian-free approximate Newton-Krylov startup strategy for RANS simulations. - Journal of Computational Physics, 397, 108741. - - [9] - Secco, N., Kenway, G. K. W., He, P., Mader, C. A., and Martins, J. R. R. A. (2021). - Efficient Mesh Generation and Deformation for Aerodynamic Shape Optimization. - AIAA Journal. - - [10] - Wu, N., Kenway, G., Mader, C. A., Jasa, J., and Martins, J. R. R. A. (2020). - pyOptSparse: A Python framework for large-scale constrained nonlinear optimization of sparse systems. - Journal of Open Source Software, 5(54), 2564. - - [11] - Hajdik, H. M., Yildirim, A., Wu, N., Brelje, B. J., Seraj, S., Mangano, M., Anibal, J. L., Jonsson, E., Adler, E. J., Mader, C. A., Kenway, G. K. W., and Martins, J. R. R. A. (2023). - pyGeo: A geometry package for multidisciplinary design optimization. - Journal of Open Source Software, 8(87), 5319. - - ## Lead - Cashen Diniz @cashend """ version = 0 @@ -275,14 +110,17 @@ class Conditions: mach: Annotated[ float, bounded(lower=0.0).category(IMPL), bounded(lower=0.1, upper=1.0).warning().category(IMPL) ] = 0.8 + """Mach number""" reynolds: Annotated[ float, bounded(lower=0.0).category(IMPL), bounded(lower=1e5, upper=1e9).warning().category(IMPL) ] = 1e6 + """Reynolds number""" area_initial: float = float("NAN") """actual initial airfoil area""" area_ratio_min: Annotated[float, bounded(lower=0.0, upper=1.2).category(THEORY)] = 0.7 """Minimum ratio the initial area is allowed to decrease to i.e minimum_area = area_initial*area_target""" cl_target: float = 0.5 + """Target lift coefficient to satisfy equality constraint""" conditions = Conditions() diff --git a/engibench/problems/beams2d/v0.py b/engibench/problems/beams2d/v0.py index d6509783..c25c7dcd 100644 --- a/engibench/problems/beams2d/v0.py +++ b/engibench/problems/beams2d/v0.py @@ -54,116 +54,10 @@ def volume_fraction_bound(design: npt.NDArray, volfrac: float) -> None: class Beams2D(Problem[npt.NDArray]): r"""Beam 2D topology optimization problem. - ## Problem Description Beams2D is a structural topology optimization (TO) problem that optimizes a 2D Messerschmitt-Bölkow-Blohm (MBB) beam under bending. The beam is symmetric about the central vertical axis, with a force applied at the top; only the right half is modeled in our case. Problems are formulated using density-based TO, drawing from an existing Python [implementation](https://github.com/arjendeetman/TopOpt-MMA-Python). - - ## Motivation - The optimization of beam cross-sections is one of a fundamental problem in engineering, aiming to - maximize the structural stiffness under some applied force. This objective is usually formulated as - minimizing the compliance, which is the inverse of stiffness. In particular, TO frames the problem as - one of optimal material distribution, defining a grid of elements for which the material densities must - be determined on a scale from 0 to 1, where 1 represents the presence of material. After applying the - beam loads and other boundary conditions, designs are typically optimized using a gradient-based - approach with the help of the finite element method (FEM). While this is one of the simplest - TO applications, it is still a computationally expensive process requiring many iterations, opening the - door for faster approximation methods such as generative inverse design. - One of the most common beam types in TO is the Messerschmitt-Bölkow-Blohm (MBB) beam, - which is supported at the bottom-right and bottom-left corners, with a downward force applied - on the top-center. Given this symmetric configuration, one half of the design may be optimized - while representing the entire structure. We implement the MBB beam in ENGIBENCH for the most - accessible comparison to previous works in this domain. - - ## Design Space - This problem simulates the right half-section of a MBB beam under bending. This half-beam is - subjected to a force at its top-left corner (corresponding to the top-center of the entire design) which - may also be shifted to the right to simulate different loading conditions. A roller support at the - bottom-right corner prevents vertical movement, and a symmetric boundary condition is enforced on - the left edge. The design space is an array of solid densities in `[0., 1.]` with a default size of - `(100, 50)` used by default, where `nelx = 100` and `nely = 50`. Internally, this is represented as - a flattened `(5000,)` array. Alternative shapes include `(50, 25)` for faster computation and `(200, 100)` - for higher-resolution results. Corresponding datasets for these three resolutions are provided. - - ## Objectives - The goal is to optimize the distribution of solid material to minimize compliance - (equivalently, maximize stiffness) while satisfying constraints on material usage and minimum feature size. - Compliance is calculated as the sum of strain energy over the structure. - - The objectives are defined and indexed as follows: - - 0. `c`: Compliance to minimize. - - ## Conditions - The following input parameters define the problem conditions: - - `volfrac`: Desired volume fraction of solid material. - - `rmin`: Minimum feature length of beam members. - - `forcedist`: Fractional distance of the downward force from the top-left (default) to the top-right corner. - - `overhang_constraint`: Boolean flag to enable a 45-degree overhang constraint for manufacturability. - - ## Simulator - Our simulation code is based on a Python adaptation of the popular 88-line topology optimization - code. It uses the more versatile density filtering approach in combination with a standard - Optimality Criteria (OC) optimization method. Two primary sensitivity matrices, one with respect - to compliance (`dc`) and the other with respect to volume fraction (`dv`), are continuously updated - and used to calculate a given design's compliance value. We have also ensured that during the - required Lagrange multiplier search within OC, the inner optimization loop terminates if the absolute - difference upper and lower bounds diminishes to a value smaller than machine precision. This - prevents the code from becoming stuck at this point, which we observed in some warm-starting - instances with noisy initial designs. - - Compliance `c` is calculated using: - ```python - c = ((Emin + xPrint**penal * (Emax - Emin)) * ce).sum() - ``` - - where `xPrint` is the current true density field, `penal` is the penalization factor (e.g., 3.0), - and `ce` is the element-wise strain energy density. - - ## Dataset - This problem offers multiple datasets for various sizes of `nelx` and `nely`. Each dataset includes - columns for the optimal design, all conditions listed above, and the corresponding objective values. - For advanced usage, we also provide a column containing the optimization history. The datasets have - been generated by sampling conditions over a structured grid for various problem sizes. - Three datasets are available on the [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab). - They correspond to resolutions of $50 \times 25$, $100 \times 50$ (default), and $200 \times 100$. - - ### v0 - - #### Fields - Each dataset contains: - - Optimized beam structures, - - The corresponding condition parameters, - - Objective values (compliance), - - Full optimization histories (for advanced use). - - #### Creation Method - Datasets were generated by uniformly sampling the condition space. The resolutions used are: - - `(50, 25)` - - `(100, 50)` - - `(200, 100)` - - A more comprehensive description of the creation method can be found in the [README](https://github.com/IDEALLab/EngiBench/tree/main/engibench/problems/beams2d). - - ## Citation - If you use this problem in your research, please cite the following paper: - ``` - @article{andreassen2011efficient, - title={Efficient topology optimization in MATLAB using 88 lines of code}, - author={Andreassen, Erik and Clausen, Anders and Schevenels, Mattias and Lazarov, Boyan S and Sigmund, Ole}, - journal={Structural and Multidisciplinary Optimization}, - volume={43}, - number={1}, - pages={1--16}, - year={2011}, - publisher={Springer} - } - ``` - - - ## Lead - Arthur Drake @arthurdrake1 """ version = 0 @@ -178,11 +72,15 @@ class Conditions: bounded(lower=0.0, upper=1.0).category(THEORY), bounded(lower=0.1, upper=0.9).warning().category(IMPL), ] = 0.35 + """Desired volume fraction of solid material""" rmin: Annotated[ float, greater_than(0.0).category(THEORY), bounded(lower=1.0, upper=10.0).category(IMPL).warning() ] = 2.0 + """Minimum feature length of beam members""" forcedist: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.0 + """Fractional distance of the downward force from the top-left (default) to the top-right corner""" overhang_constraint: bool = False + """Boolean flag to enable a 45-degree overhang constraint for manufacturability""" @dataclass class SimulateConfig(Conditions): diff --git a/engibench/problems/beams2d/v1.py b/engibench/problems/beams2d/v1.py index 089d5657..8965d96c 100644 --- a/engibench/problems/beams2d/v1.py +++ b/engibench/problems/beams2d/v1.py @@ -26,13 +26,12 @@ class Beams2D(Beams2D_v0): r"""Beam 2D topology optimization problem - Version 1 (v1). - ### v1 + ## v1 This version augments v0 by fixing a minor detail in the v0 warm-start optimization process. Specifically, when warm-starting from a provided design, a small epsilon value is added to avoid zero-density values that could lead to gradient issues. The datasets themselves remain unchanged. All other behavior is identical to v0. - See v0.py for full baseline documentation. """ version = 1 diff --git a/engibench/problems/heatconduction2d/v0.py b/engibench/problems/heatconduction2d/v0.py index d0860f6d..d69c4457 100644 --- a/engibench/problems/heatconduction2d/v0.py +++ b/engibench/problems/heatconduction2d/v0.py @@ -40,50 +40,12 @@ def volume_fraction_bound(design: npt.NDArray, volume: float) -> None: class HeatConduction2D(Problem[npt.NDArray]): r"""HeatConduction 2D topology optimization problem. - ## Problem Description This problem simulates the performance of a Topology optimisation of heat conduction problems governed by the Poisson equation (https://www.dolfin-adjoint.org/en/stable/documentation/poisson-topology/poisson-topology.html) ## Motivation Heat conduction problems serve as fundamental benchmarks for the development and evaluation of design optimization methods, with applications ranging from thermal management in electronic devices to insulation systems and heat exchangers in industrial applications. As thermal management has become critical in fields such as aerospace, automotive, and consumer electronics, both industry and academia have shown growing interest in advanced thermal design systems. In response to this demand, topology optimization has become popular as a powerful approach for improving heat dissipation while minimizing material usage. In addition, the development of additive manufacturing technologies has made the complex geometries produced by topology optimization more feasible to fabricate in real-world applications. - - ## Design space - These problems are a specific subset of topology optimization problems aimed at minimizing thermal compliance within a unit square (2D) subject to: a constraint on the volume of highly conductive material used, and given boundary conditions, - particularly the location of the adiabatic region. The adiabatic region refers to a symmetric length on the bottom side of the 2D problem space . The design space for the 2D problem consists of a 2D array representing solid densities, - which is parametrized by resolution, that is, 'DesignSpace = [0,1] By default, a $101 \times 101$ space is used for the 2D problem. - - ## Objectives - The objective is defined and indexed as follows: - - 0. `c`: Thermal compliance coefficient to minimize. - - ## Conditions - The conditions are defined by the following parameters: - - `volume`: the volume limits on the material distributions - - `length`: The length of the adiabatic region on the bottom side of the design domain. - - ## Simulator - The simulator is a docker container with the dolfin-adjoint software that computes the thermal compliance of the design. - We convert use intermediary files to convert from and to the simulator that is run from a Docker image. - - ## Dataset - The dataset has been generated the dolfin-adjoint software. It is hosted on the [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/heat_conduction_2d_v0). - - ### v0 - - #### Fields - The dataset only contains conditions and optimal designs (no objective). - - #### Creation Method - The creation method for the dataset is specified in the reference paper. - - ## References - If you use this problem in your research, please cite the following paper: - Milad Habibi, Jun Wang, and Mark Fuge, "When Is it Actually Worth Learning Inverse Design?" in IDETC 2023. doi: https://doi.org/10.1115/DETC2023-116678 - - ## Lead - Milad Habibi @MIladHB """ version = 0 @@ -98,9 +60,9 @@ class Conditions: bounded(lower=0.0, upper=1.0).category(THEORY), bounded(lower=0.3, upper=0.6).warning().category(IMPL), ] = 0.5 - """Volume constraint""" + """Volume limits on the material distributions""" length: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.5 - """Length constraint""" + """Length of the adiabatic region on the bottom side of the design domain""" @dataclass class Config(Conditions): diff --git a/engibench/problems/heatconduction3d/v0.py b/engibench/problems/heatconduction3d/v0.py index bbababef..5bb823a1 100644 --- a/engibench/problems/heatconduction3d/v0.py +++ b/engibench/problems/heatconduction3d/v0.py @@ -40,51 +40,8 @@ def volume_fraction_bound(design: npt.NDArray, volume: float) -> None: class HeatConduction3D(Problem[npt.NDArray]): r"""HeatConduction 3D topology optimization problem. - ## Problem Description + Problem Description This problem simulates the performance of a Topology optimisation of heat conduction problems governed by the Poisson equation (https://github.com/dolfin-adjoint/pyadjoint/blob/master/examples/poisson-topology/poisson-topology.py) - - ## Motivation - Heat conduction problems serve as fundamental benchmarks for the development and evaluation of design optimization methods, with applications ranging from thermal management in electronic devices to insulation systems and - heat exchangers in industrial applications. As thermal management has become critical in fields such as aerospace, automotive, and consumer electronics, both industry and academia have shown growing interest in advanced thermal - design systems. In response to this demand, topology optimization has become popular as a powerful approach for improving heat dissipation while minimizing material usage. In addition, the development of additive manufacturing - technologies has made the complex geometries produced by topology optimization more feasible to fabricate in real-world applications. - - ## Design space - These problems are a specific subset of topology optimization problems aimed at minimizing thermal compliance within a unit cube (3D), subject to: a constraint on the volume of highly conductive material used, and given boundary conditions, - particularly the location of the adiabatic region. The adiabatic region refers to a prescribed symmetric area on the bottom surface of the 3D problem space. The 3D design space is represented as a 3D tensor of densities. By default, - a $51 \times 51 \times 51$ space is used for the 3D problem. - - ## Objectives - The objective is defined and indexed as follows: - - 0. `c`: Thermal compliance coefficient to minimize. - - ## Conditions - The conditions are defined by the following parameters: - - `volume`: the volume limits on the material distributions - - `area`: The area of the adiabatic region on the bottom side of the design domain. - - ## Simulator - The simulator is a docker container with the dolfin-adjoint software that computes the thermal compliance of the design. - We convert use intermediary files to convert from and to the simulator that is run from a Docker image. - - ## Dataset - The dataset has been generated the dolfin-adjoint software. It is hosted on the [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/heat_conduction_3d_v0). - - ### v0 - - #### Fields - The dataset only contains conditions and optimal designs (no objective). - - #### Creation Method - The creation method for the dataset is specified in the reference paper. - - ## References - If you use this problem in your research, please cite the following paper: - Habibi, Milad, Shai Bernard, Jun Wang, and Mark Fuge, "Mean squared error may lead you astray when optimizing your inverse design methods" in JMD 2025. doi: https://doi.org/10.1115/1.4066102 - - ## Lead - Milad Habibi @MIladHB """ version = 0 @@ -99,9 +56,9 @@ class Conditions: bounded(lower=0.0, upper=1.0).category(THEORY), bounded(lower=0.3, upper=0.6).warning().category(IMPL), ] = 0.3 - """Volume constraint""" + """Volume limits on the material distributions""" area: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.5 - """Area constraint""" + """Area of the adiabatic region on the bottom side of the design domain""" @dataclass class Config(Conditions): diff --git a/engibench/problems/photonics2d/v0.py b/engibench/problems/photonics2d/v0.py index 8ed3196f..79607004 100644 --- a/engibench/problems/photonics2d/v0.py +++ b/engibench/problems/photonics2d/v0.py @@ -53,162 +53,10 @@ class Photonics2D(Problem[npt.NDArray]): r"""Photonic Inverse Design 2D Problem (Wavelength Demultiplexer). - ## Problem Description Optimize a 2D material distribution (`rho`) to function as a wavelength demultiplexer, routing wave with `lambda1` to output 1 and `lambda2` to output 2. The design variables represent material density that is converted to permittivity using filtering and projection. - - ### Motivation - The optimization of photonic circuits in general, and multiplexers in particular, was - one of the initial and most widely studied problems in the inverse design of - electromagnetic/optical devices. In part, this is because multiplexer - devices have several interesting properties that make them more difficult to create generative - models of, compared to other problems in EngiBench. This includes the fact that, due to the wave - properties of the physical phenomena there are usually multiple solutions with equivalent or - similar performance, which results from shifting or inverting the phase profile of the - electromagnetic wave. This adds complexity to the generative model in that the solution may not - have a single unique global minimum. Another motivating factor for including this problem is the - complexity of the structures/designs themselves: unlike in structural or thermal compliance - problems, which lead to connected structures, the photonics solutions often involve several - disconnected elements whose relative position and spacing is governed by the specific wavelengths - it needs to demultiplex. This is a difficult prediction and generation task, compared to, e.g., - generating a connected beam structure. Thus, it acts as a good counterpoint to add to the library - and provides a mechanism to benchmark generative algorithms that can perform well on both - connected and disconnected design topologies. - - ## Design space - This problem simulates a wavelength demultiplexer where the optimized device will direct an - electromagnetic wave to one of two possible output ports depending on the wavelength/frequency of - the incoming wave. Specifically, the demultiplexer targets two specific wavelengths (referred to - $\lambda_1$ and $\lambda_2$ in the library), and the performance of the device is how well it can - bend or direct the energy toward two specific locations in the device, as measured by how much of - the electric field of each wavelength overlaps with the desired output port locations. The design - space consists of a 2D array representing the presence of either a high or low permittivity, - parameterized by `nelx` and `nely`, i.e., $\text{design_space} = [0,1]^{\text{nelx}\times \text{nely}}$. - By default, the library uses a $120 \times 120$ space, however, this can be modified to non-square - design spaces by the user. Specifically, we use a 2D tensor `rho` (num_elems_x, num_elems_y) with - values in [0, 1], representing material density. This is stored as `design_space` (gymnasium.spaces.Box). - - ## Objectives - The main objective is to maximize the overlap of the electric field of the simulated wavelength at - the target output location, with an optional penalty for the amount of material used (this penalty - weight is set to a small default value ($1e^{-2}$) for consistency, but can be altered for advanced - usage): - 0. `total_overlap`: Objective to maximize, defined as - `overlap1 * overlap2`. Higher is better. This corresponds to - the overlap in the target electrical fields with the desired demultiplexing locations. - Note that bot `simulate` and `optimize` subtract a small material penalty - (`total_overlap - penalty`) to avoid multiple equivalent local optima, but this penalty - is small relative to the overlap objective. - - ## Conditions - These are designed as user-configurable parameters that alter the problem definition. The - conditions consist of the two input wavelengths to be demultiplexed -- - $\lambda_1$ and $\lambda_2$, as well as a desired `blur_radius` ($r_{blur}$) parameter, - which blurs (using a circular convolution) the pixelized design field for a chosen number of - integer pixels\textemdash this blurring essentially creates a penalty on the minimum feature size - of the design. The size of the device -- expressed as `nelx` and `nely` -- is also adjustable, - and could be viewed as a possible condition for multi-resolution problems, but in practice, as with - `Beams2D`, this is built into the problem definition since it produces a different dataset. - - Default problem parameters that can be overridden via the `config` dict: - - `lambda1`: The first input wavelength in μm (default: 1.5 μm). - - `lambda2`: The first input wavelength in μm (default: 1.3 μm). - - `blur_radius`: Radius for the density blurring filter (default: 2). - Higher values correspond to larger elements, which could - possibly be more manufacturable. - - `num_elems_x`: Number of grid cells in x (default: 120). - - `num_elems_y`: Number of grid cells in y (default: 120). - - In practice, for the dataset loading, we will keep `num_elems_x` and `num_elems_y`to set - values for each dataset, such that different resolutions correspond to different - independent datasets. - - ## Optimization Parameters - Note: These are advanced parameters that alter the optimization process -- - we do not recommend changing these if you are only using the library for benchmarking, - as it could make results less reproducible across papers using this problem.) - - `num_optimization_steps`: Total number of optimization steps (default: 300). - - `step_size`: Adam optimizer step size (default: 1e-1). - - `penalty_weight`: Weight for the L2 penalty term (default: 1e-2). Larger values reduce - unnecessary material, but may lead to worse performance if too large. - - `eta`: Projection center parameter (default: 0.5). There is little reason to change this. - - `N_proj`: Number of projection applications (default: 1). Increasing this can help make - the design more binary. - - `N_blur`: Number of blur applications (default: 1). Increasing this smooths the design more. - - `initial_beta`: Initial beta for the optimization continuation scheme (default: 1.0). - - `save_frame_interval`: Interval for saving intermediate design frames during optimization. - If > 0, saves a frame every `save_frame_interval` iterations - to the `opt_frames/` directory. Default is 0 (disabled). - - ## Internal Constants - Note: These are not typically changed by users, but provided here for technical reference - - `dl`: Spatial resolution (meters) (default: 40e-9). - - `Npml`: Number of PML cells (default: 20). - - `epsr_min`: Minimum relative permittivity (default: 1.0). - - `epsr_max`: Maximum relative permittivity (default: 12.0). - - `space_slice`: Extra space for source/probe slices (pixels) (default: 8). - - ## Simulator - The simulation uses the `ceviche` library's Finite Difference Frequency Domain (FDFD) - solver (`fdfd_ez`). Optimization uses `ceviche.optimizers.adam_optimize` with - gradients computed via automatic differentiation (`autograd`). - The simulation code uses the \texttt{ceviche} library and specifically, - the wave demultiplexer demonstration case provided by the - [library authors](https://github.com/fancompute/workshop-invdesign/blob/master/04_Invdes_wdm_scheduling.ipynb} - based on their related publication, which uses a similar formalism to an earlier demultiplexer - paper by Piggott (2015). The optimization method is first-order and uses the Adam optimizer. Beyond - the baseline implementation already available via `ceviche`, we implemented a polynomial $\beta$ continuation - scheme that performed more reliably than the step-wise continuation scheme used in the original implementation, - and EngiBench also possesses the ability to change the starting and ending continuation values, for future - research cases where one wishes to estimate or optimize the continuation schedule themselves. Other than these - changes, the implementation of this problem is as consistent as possible with that of the original `ceviche` library. - - ## Dataset - This problem offers a single datasets of `nelx`=120 and `nely`=120, although various sizes of `nelx` and `nely` - could be generated from the library if desired. The dataset includes columns for the optimal design, all conditions - listed above, and the corresponding objective value history as the optimizer optimized toward the optimal design - provided in the dataset. The dataset was generated by sampling by sampling at random $\lambda_1$, $\lambda_2$ and - $r_{blur}$ over the conditions mentioned above. The dataset is available on the - [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/photonics_2d_120_120_v0). - - ### v0 - - #### Fields - Each dataset contains: - - `lambda1`: The first input wavelength in μm. - - `lambda2`: The second input wavelength in μm. - - `blur_radius`: Radius for the density blurring filter (pixels). - - `optimal_design`: The optimal design density array (shape num_elems_x, num_elems_y). - - `optimization_history`: A list of objective values from the optimization process (field overlap minus penalty, where higher is better) -- This is for advanced use. - - #### Creation Method - To generate a dataset for training, we generate (randomly, uniformly) swept over the following parameters: - - $\lambda_1 \in [0.5\mu m, 1.25\mu m]$ = `lambda1` = `rng.uniform(low=0.5, high=1.25, size=20)` -- This corresponds roughly to a portion of the visible spectrum up to near-infrared. - - $\lambda_2 \in [0.75\mu m, 1.5\mu m]$ = `lambda2` = `rng.uniform(low=0.75, high=1.5, size=20)` -- This corresponds roughly to a portion of the visible spectrum up to near-infrared. - - $r_{blur}$ = `blur_radius` = `range(0, 5)` - - ## Citation - This problem is directly refactored from the Ceviche Library: - https://github.com/fancompute/ceviche - and if you use this problem your experiments, you can use the citation below - provided by the original library authors: - ``` - @article{hughes2019forward, - title={Forward-Mode Differentiation of Maxwell's Equations}, - author={Hughes, Tyler W and Williamson, Ian AD and Minkov, Momchil and Fan, Shanhui}, - journal={ACS Photonics}, - volume={6}, - number={11}, - pages={3010--3016}, - year={2019}, - publisher={ACS Publications} - } - ``` - - ## Lead - Mark Fuge @markfuge """ version = 0 @@ -258,7 +106,8 @@ class Conditions: blur_radius: Annotated[ int, bounded(lower=0).category(THEORY | IMPL), bounded(lower=0, upper=5).warning().category(IMPL) ] = 2 - """Radius for the density blurring filter (pixels)""" + """Radius for the density blurring filter (pixels). + Higher values correspond to larger elements, which could possibly be more manufacturable.""" @dataclass class Config(Conditions): diff --git a/engibench/problems/power_electronics/v0.py b/engibench/problems/power_electronics/v0.py index 2f2a4c9c..62c31b2c 100644 --- a/engibench/problems/power_electronics/v0.py +++ b/engibench/problems/power_electronics/v0.py @@ -31,90 +31,10 @@ class PowerElectronics(Problem[npt.NDArray]): This problem requires `ngspice` to be installed. See the simulator section for more details. ``` - ## Problem Description This problem simulates a power converter circuit which has a fixed circuit topology. There are 5 switches, 4 diodes, 3 inductors and 6 capacitors. The circuit topology is fixed. It is defined in the netlist file `5_4_3_6_10-dcdc_converter_1.net`. By changing circuit parameters such as capacitance, we rewrite the netlist file and use ngSpice to simulate the circuits to get the performance metrics, which are defined as the objectives of this problem. You can use this problem to train your regression model. You can also try to find the optimal circuit parameters that minimize objectives. - - ## Motivation - Optimizing circuit parameters is a critical aspect of circuit design but remains challenging, particularly for power converter circuits that contain diodes and switches, - which introduce significant nonlinearity and discontinuity. These characteristics make key objectives such as *DcGain* and *Voltage Ripple* highly sensitive to even small parameter variations. - - Because the circuit simulator NgSpice operates as a black box and is non-differentiable, gradient-based optimization methods are not suitable. - Bayesian optimization is commonly employed for parameter tuning, while surrogate models offer a promising alternative. - Even under the constraint of a fixed topology, optimizing circuit parameters to minimize the objectives remains a difficult problem for surrogate models. - - NgSpice applies transient analysis by formulating the system as a set of differential equations based on Kirchhoff's laws. - These equations are discretized using numerical integration methods such as the Backward Euler or trapezoidal rule and solved iteratively at each time step to compute performance metrics. - To ensure stable simulations, a specific on-off switching pattern is chosen for the circuit. - Despite this simplification, determining the optimal parameter values remains highly challenging. - - ## Design Space - The design space for this problem is represented as a 10-dimensional bounded box, where each dimension corresponds to a specific circuit parameter. These parameters include values for capacitors, inductors, and a shared duty cycle for all switches. Each design can be expressed as a vector **x** of the form: - - $$ - x = \begin{bmatrix} C_1,\dots,C_6,L_1,L_2,L_3,T_1 \end{bmatrix}^{\top} \in \mathcal{X}, - \quad - \mathcal{X} = [1\text{e}{-6}, 2\text{e}{-5}]^6 \times [1\text{e}{-6}, 1\text{e}{-3}]^3 \times [0.1, 0.9] - $$ - - Here, $C_1,\dots,C_6$ are the capacitance values (in Farads), $L_1,L_2,L_3$ are the inductance values (in Henries), and $T_1$ is the duty cycle shared across all 5 switches. The duty cycle $T_1$ denotes the fraction of time during which the switches are in the “on” state and governs a periodic on-off pattern repeated at high frequency throughout the simulation. - - ## Objectives - The simulation outputs two scalar values: *DcGain* and *Voltage Ripple*. The former represents the ratio of load to input voltage and should ideally approximate a predefined constant, such as $0.25$, as closely as possible. Meanwhile, the latter quantifies the voltage fluctuation at the load. - - **DcGain objective:** - - $$ - \min_{\mathbf{x} \in \mathcal{X}} \; \bigg|\frac{\overline{V_{load}(t)}}{V_{source}} - 0.25\bigg| - = \bigg|\frac{1}{V_{source}} \cdot \frac{1}{T} \sum_{i=1}^{N-1} \frac{V_{load}(t_{i+1}) + V_{load}(t_i)}{2} \cdot (t_{i+1} - t_i) - 0.25\bigg| - $$ - - where $\overline{V_{load}(t)}$ is the average load voltage, $V_{source} = 1000$ volts, and $T = t_N - t_1$ is the simulation duration. - - **Voltage Ripple objective:** - - $$ - \min_{\mathbf{x} \in \mathcal{X}} \; \text{Voltage Ripple} - = \frac{V_{pp}(t)}{\overline{V_{load}(t)}} - = \frac{\max_{i \in [1, N]} V_{load}(t_i) - \min_{i \in [1, N]} V_{load}(t_i)}{\overline{V_{load}(t)}} - $$ - - where $V_{pp}$ is the peak-to-peak load voltage calculated during transient analysis. - - ## Conditions - This problem does not include environmental or operational conditions as part of its input specification. Unlike other domains where the simulation setup may vary based on conditions (e.g., load configurations or external temperatures), the circuit is simulated under fixed source voltage and switching behavior. As a result, the design optimization task focuses solely on tuning internal circuit parameters, with no external conditions to vary. More complex variants of this problem — involving multiple topologies or variable source voltages — may be considered in future releases. - - ## Simulator - The simulator is ngSpice circuit simulator. You can download it based on your operating system: - - Windows: [https://sourceforge.net/projects/ngspice/files/ng-spice-rework/45.2/](https://sourceforge.net/projects/ngspice/files/ng-spice-rework/45.2/) - - MacOS: `brew install ngspice` - - Linux: `sudo apt-get install ngspice` - - ## Dataset - The dataset linked to this problem is hosted on the [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/power_electronics). - - ### v0 - - #### Fields - The dataset contains 3 fields: - - `initial_design`: The 20-dimensional design variable defined above. - - `DcGain`: The ratio of load vs. input voltage. - - `Voltage_Ripple`: The fluctuation of voltage on the load `R0`. - - #### Creation Method - We created this dataset in 3 parts. All the 3 parts are simulated with {`GS0_L1`, `GS1_L1`, `GS2_L1`, `GS3_L1`, `GS4_L1`} = {1, 0, 0, 1, 1} and {`GS0_L2`, `GS1_L2`, `GS2_L2`, `GS3_L2`, `GS4_L2`} = {1, 0, 1, 1, 0}. - Here are the 3 parts: - 1. 6 capacitors and 3 inductors only take their min and max values. `T1` ranges {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9}. There are 2^6 * 2^3 * 9 = 4608 samples. - 2. Random sample 4608 points in the 6 + 3 + 1 = 10 dimensional space. Min and max values in each dimension will not be sampled. - 3. Latin hypercube sample 4608 points in the 6 + 3 + 1 = 10 dimensional space. Each dimension is split into 10 intervals. Min and max values in each dimension will not be sampled. - - ## References - If you use this problem in your research, please cite the following paper: - - ## Lead - Xuliang Dong @ liangXD523 """ version = 0 diff --git a/engibench/problems/thermoelastic2d/v0.py b/engibench/problems/thermoelastic2d/v0.py index bdbe7cd3..de42096d 100644 --- a/engibench/problems/thermoelastic2d/v0.py +++ b/engibench/problems/thermoelastic2d/v0.py @@ -34,65 +34,7 @@ class ThermoElastic2D(Problem[npt.NDArray]): r"""Truss 2D integer optimization problem. - ## Problem Description This is 2D topology optimization problem for minimizing weakly coupled thermo-elastic compliance subject to boundary conditions and a volume fraction constraint. - - ## Motivation - As articulated in their respective sections, both the Beams2D and HeatConduction2D problems found in the EngiBench library are fundamental engineering design problems that have historically served as benchmarks for the development and testing of optimization methods. - While their relevance is supported by needs in real engineering design scenarios (aerospace, automotive, consumer electronics, etc...), their mono-domain nature ignores the reality that coupling between domains exists, and should be accounted for in scenarios where performance in one domain significantly impacts performance in another. - To address this distinction, a multi-physics topology optimization problem is developed that captures the coupling between structural and thermal domains. - - ## Design space - This multi-physics topology optimization problem is governed by linear elasticity and steady-state heat conduction with a one-way coupling from the thermal domain to the elastic domain. - The problem is defined over a square 2D domain, where load elements and support elements are placed along the boundary to define a unique elastic condition. - Similarly, heatsink elements are placed along the boundary to define a unique thermal condition. - The design space is then defined by a 2D array representing density values (parameterized by DesignSpace = [0,1]^{nelx x nely}, where nelx and nely denote the x and y dimensions). - - ## Objectives - The objective of this problem is to minimize total compliance C under a volume fraction constraint V by placing a thermally conductive material. - Total compliance is defined as the sum of thermal compliance and structural compliance. - - ## Conditions - Problem conditions are defined by creating a python dict with the following info: - - `fixed_elements`: Encodes a binary NxN matrix of the structurally fixed elements in the domain. - - `force_elements_x`: Encodes a binary NxN matrix specifying elements that have a structural load in the x-direction. - - `force_elements_y`: Encodes a binary NxN matrix specifying elements that have a structural load in the y-direction. - - `heatsink_elements`: Encodes a binary NxN matrix specifying elements that have a heat sink. - - `volfrac`: Encodes the target volume fraction for the volume fraction constraint. - - `rmin`: Encodes the filter size used in the optimization routine. - - `weight`: Allows one to control which objective is optimized for. 1.0 Is pure structural optimization, while 0.0 is pure thermal optimization. - - ## Simulator - The simulation code is based on a Python adaptation of the popular 88-line topology optimization code, modified to handle the thermal domain in addition to thermal-elastic coupling. - Optimization is conducted by reformulating the integer optimization problem as a continuous one (leveraging a SIMP approach), where a density filtering approach is used to prevent checkerboard-like artifacts. - The optimization process itself operates by calculating the sensitivities of the design variables with respect to total compliance (done efficiently using the Adjoint method), calculating the sensitivities of the design variables with respect to the constraint value, and then updating the design variables by solving a convex-linear subproblem and taking a small step (using the method of moving asymptotes). - The optimization loop terminates when either an upper bound of the number of iterations has been reached or if the magnitude of the gradient update is below some threshold. - - ## Dataset - The dataset linked to this problem is on huggingface [Hugging Face Datasets Hub](https://huggingface.co/datasets/IDEALLab/thermoelastic_2d_v0). - This dataset contains a set of 1000 optimized thermoelastic designs in a 64x64 domain, where each design is optimized for a unique set of conditions. - Each datapoint's conditions are randomly generated by arbitrarily placing: a single loaded element along the bottom boundary, two fixed elements (fixed in both the x and y direction) along the left and top boundary, and heatsink elements along the right boundary. - Furthermore, values for the volume fraction constraint are randomly selected in the range $[0.2, 0.5]$. - - Relevant datapoint fields include: - - `optimal_design`: An optimized design for the set of boundary conditions - - `strain`: The strain field for the initial uniform design - - `vm_stress`: The von Mises stress field for the initial uniform design - - `fixed_elements`: Encodes a binary NxN matrix of the structurally fixed elements in the domain. - - `force_elements_x`: Encodes a binary NxN matrix specifying elements that have a structural load in the x-direction. - - `force_elements_y`: Encodes a binary NxN matrix specifying elements that have a structural load in the y-direction. - - `heatsink_elements`: Encodes a binary NxN matrix specifying elements that have a heat sink. - - `volume_fraction`: The volume fraction value of the optimized design - - `structural_compliance`: The structural compliance of the optimized design - - `thermal_compliance`: The thermal compliance of the optimized design - - `nelx`: The number of elements in the x-direction - - `nely`: The number of elements in the y-direction - - `volfrac`: The volume fraction target of the optimized design - - `rmin`: The filter size used in the optimization routine - - `weight`: The domain weighting used in the optimization routine - - ## Lead - Gabriel Apaza @gapaza """ version = 0 @@ -109,21 +51,27 @@ class Conditions: fixed_elements: Annotated[npt.NDArray[np.int64], bounded(lower=0.0, upper=1.0).category(THEORY)] = field( default_factory=lambda: FIXED_ELEMENTS ) + """Binary NxN matrix of the structurally fixed elements in the domain""" force_elements_x: Annotated[npt.NDArray[np.int64], bounded(lower=0.0, upper=1.0).category(THEORY)] = field( default_factory=lambda: FORCE_ELEMENTS_X ) + """Binary NxN matrix specifying elements that have a structural load in the x-direction""" force_elements_y: Annotated[npt.NDArray[np.int64], bounded(lower=0.0, upper=1.0).category(THEORY)] = field( default_factory=lambda: FORCE_ELEMENTS_Y ) + """Binary NxN matrix specifying elements that have a structural load in the y-direction""" heatsink_elements: Annotated[npt.NDArray[np.int64], bounded(lower=0.0, upper=1.0).category(THEORY)] = field( default_factory=lambda: HEATSINK_ELEMENTS ) + """Binary NxN matrix specifying elements that have a heat sink""" volfrac: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.3 + """Target volume fraction for the volume fraction constraint""" rmin: Annotated[ float, bounded(lower=1.0).category(THEORY), bounded(lower=0.0, upper=3.0).warning().category(IMPL) ] = 1.1 + """Filter size used in the optimization routine""" weight: Annotated[float, bounded(lower=0.0, upper=1.0).category(THEORY)] = 0.5 - """1.0 for pure structural, 0.0 for pure thermal""" + """Control which objective is optimized for. 1.0 is pure structural optimization, while 0.0 is pure thermal optimization""" conditions = Conditions() design_space = spaces.Box(low=0.0, high=1.0, shape=(NELX, NELY), dtype=np.float32)