Efficient neural topology optimization via active learning for enhancing turbulent mass transfer in fluid channels

We first investigate the effects of topology and inlet velocity on turbulent mass transfer. The concentration and pressure distributions for the smooth channel and channels with randomly generated topological structures (Figs. 1b and c) are numerically simulated based on the mechanism model presented in Section 4.1. Comparisons between the results of the three fluid channels (Fig. 1c) show that modifying the flow channel by adding solid blocks generally results in a more uniform concentration distribution (smaller $\sigma_{C}^{2}$) but with the cost of increase of the system pressure drop (larger $\Delta P$) compared to the smooth channel. The Structure ii achieves smaller values of both $\Delta P$ and $\sigma_{C}^{2}$ compared to that of Structure i, and therefore Structure ii is more favorable for the mass transfer process. This demonstrates that different channel structures behave differently in terms of the two criteria, and optimization is necessary to find the best topological structure of the channel that minimizes the overall objective function.

The inlet velocity of the system is also a crucial factor influencing mass transfer efficiency. The system pressure drop $\Delta P$ increases with the increase of inlet velocity (Fig. 1d), whereas the $\sigma_{C}^{2}$ decreases first and then increases with the inlet velocity. The concentration distribution is the most uniform at $v=0.5\,\mathrm{m/s}$. This suggests that increasing the inlet velocity improves both the convective and fluctuating mass transfer in the turbulent flow system. At the same time, the flow rate of the solute entering the system also increases with the rise of inlet velocity $v$. Therefore, the high-concentration region becomes larger at $v=0.9\,\mathrm{m/s}$ compared to $v=0.1\,\mathrm{m/s}$ , leading to an increase in the concentration variance.

Since the governing equations based on the mechanism model for the turbulent mass transfer process are complex, traditional methods cannot efficiently optimize the channel structure, especially for the investigations on various inlet velocities. In the present study, we develop a computational framework to explore the optimal topology under different inlet velocities by integrating pre-trained neural operators with a neural topology (Fig. 2). To address the challenge of the large search space of the possible topological structures, we propose an active data-augmentation method to enhance the effectiveness of the neural operator-based optimization algorithm. The method can be described by the following three steps.

After introducing the methods, we present the TO results under different inlet velocities. We first compare the objective values between the cases in the training dataset and the TO results in Fig. 3a. Specifically, we show the optimization results under three weights: $\omega=0.1$, 1, and 10. These weights can be interpreted as optimization processes dominated by enhancing mass transfer performance, balanced optimization of both objective functions, and optimization driven by minimizing external energy consumption, respectively. The dashed lines in Fig. 3a represent a possible Pareto front, which is obtained by linearly interpolating the optimization results under three different weights. When $\omega=0.1$, the optimization algorithm achieves a channel structure with the smallest $\sigma_{C}^{2}$ compared to the examples in the training set, especially for lower inlet velocity. In addition, the system pressure drop increases significantly as the inlet velocity increases, which is in agreement with the simulation results in Fig. 1d. As a result, reducing $\sigma_{C}^{2}$ no longer dominates absolutely when $v=0.9\,\mathrm{m/s}$ , and the optimal topological structure is a trade-off between the two objectives.

As $\omega$ and $v$ increase, the optimization process gradually shifts from being predominantly driven by reducing $\sigma_{C}^{2}$ to being predominantly driven by reducing $\Delta P$. Consequently, the optimal points move towards the lower right corner in Fig. 3a. When $v=0.9\,\mathrm{m/s}$ , the $\Delta P$ corresponding to the optimal topology with $\omega=1$ and 10 are close to $\Delta P$ of the smooth channel. These results indicate the effectiveness of our proposed TO method, which can output a topology with smaller $\Delta P$ when the TO objective is dominated by the system pressure drop.

Next, we compare the optimized topological structures for $v=0.1\,\mathrm{m/s}$ before and after data augmentation (Fig. 3b). When $\omega=1$, the optimal topological structures generated by the TO algorithm before and after data augmentation are quite similar (Fig. 3b, second row). However, in the case of $\omega=0.1$, the optimal topological structure obtained after data augmentation shows a noticeable difference (Fig. 3b, first row). More optimization results (Fig. S5) show that for most of the examples the proposed method is successful in getting a better result after one active data augmentation, and the second data augmentation can further fine-tune the optimization results.

As discussed above, we obtained good TO results even in the first round without data augmentation. This good performance is attributed to the good prediction accuracy and generalization ability of the neural operators trained only with the initial small dataset. The neural operators can accurately predict process objectives under different inlet velocities or topological structures. In our training dataset (Appendix S1), we use $v\in\left\{0.1,0.2,\dots 0.9\right\}\,\mathrm{m/s}$. To quantify the generalization ability, we design three datasets for testing, including a standard testing dataset with $v\in\left\{0.1,0.2,\dots 0.9\right\}\,\mathrm{m/s}$, an interpolation dataset with $v\in[0.1,0.9]\,\mathrm{m/s}$, and a difficult extrapolation dataset with $v\in\left\{0.05,0.95\right\}\,\mathrm{m/s}$. The prediction errors of the neural operators for three test datasets (Table 1) are less than 1%, 2% and 5%, respectively. Hence, the trained $\mathcal{G}_{P}$ and $\mathcal{G}_{C}$ neural operators can not only predict the new topological structure, but also has good prediction accuracy for the cases whose inlet velocities are outside the training set. Moreover, although the neural operators are not directly trained with the objective values, the objective values calculated using the predicted physical fields are in good agreement with the results from the mechanism model. We present more details on neural operator validation in Appendix S2.

Moreover, we illustrate how data augmentation, particularly for the $\mathcal{G}_{P}$ training dataset, allows the TO framework to achieve better topology. The pressure difference $\Delta P$ for the data points in the initial $\mathcal{G}_{P}$ training dataset mostly ranges between 0.004 and 0.005 (Fig. 3a), while in the TO results, it is approximately 0.003 for $\omega=0.1$, 0.002 for $\omega=1$, and 0.001 for $\omega=10$. Excluding the smooth channel, the minimum value of $\Delta P$ in the initial training dataset is around 0.003. Therefore, the training dataset of the $\mathcal{G}_{P}$ has fewer data points near the optimal point, leading to lower prediction accuracy for the optimized structure. In contrast, for the $\mathcal{G}_{C}$ neural operator (Fig. 3c), where the distribution of $\sigma_{C}^{2}$ in the training dataset is more uniform, the prediction of neural operator for optimized structure is more accurate. After adding more data points to the training set based on the TO results, the number of data points with smaller objective values increases (Fig. 3c), which improves the prediction accuracy of neural operators for optimized structure.

We make a comparison of the computational cost between our proposed method and the traditional method. Although our method requires solving the computational fluid dynamics (CFD) and computational mass transfer (CMT) model equations repeatedly for various conditions to establish the training dataset, each simulation is independent and can be performed in parallel. Based on the trained neural operator, the optimal structure of 9 inlet velocities can be obtained by the optimization algorithm. When the $\omega$ in the objective function changes, it is only necessary to repeat Step 3 to get a new optimization result. In contrast, structure optimization based on mechanism models requires solving multiple CMT and Lagrange equations iteratively. For example, Jia et al. [19] employed the MMA optimization algorithm and used 31 iterations, in each of which the mechanism models equations must be solved, to obtain the optimal channel structure with fixed values of $v$ and $\omega$.

For numerical solutions with the mechanism model, we utilizes parallel computing with 48 threads on two Intel Xeon E5-2687W v4 CPUs. The neural network training is performed using an NVIDIA GeForce RTX 3090 Ti GPU. The computational cost is summarized in Table 2, which shows that the proposed method has an obvious computational advantage in optimizing multiple structures under different inlet velocities. Furthermore, as we showed in Table 1, the trained neural operator has sufficient prediction accuracy for the physical fields and objectives corresponding to any $v\in[0.1,0.9]\,\mathrm{m/s}$. Hence, if we aim to perform TO for more inlet velocities, the speedup will be even more significant.

To validate our ML results, we conduct turbulent mass transfer experiments of the smooth and optimized flow channels under $v=0.1$ m/s and $\omega=0.1$. Fig. 4a shows the diagram of the experiment setup. The simulated 2D channel was expanded into a 3D channel with a square cross-section. Three sides of the 3D channel were printed using white material, while the remaining side was fitted with a transparent window to capture the mass transfer behavior of water-methylene blue within the channel. The dimensions of the transparent window match those of the computational domain in Fig. 1a. A centrifugal pump was used to achieve continuous inlet feed, with the inlet velocity controlled by a rotor flowmeter and a bypass control valve. During the experiment, a high-speed camera continuously captured the fluid mixing conditions within the channel, and the images were processed into grayscale on a computer. By averaging the grayscale values across all frames in the animation, the time-averaged grayscale images were obtained. The time-averaged concentration distribution (Fig. 4b) was then obtained using the time-averaged grayscale image and the calibration curve. The original experimental data and calibration curves are presented in Appendix S5.

In order to quantitatively verify the effectiveness of the optimized channel in enhancing mass transfer, the variance of concentration distributions ($\sigma_{C}^{2}$) for both the smooth and the optimized fluid channels are compared in Fig. 4c. Experimental results indicate that the mass transfer performance of the optimized fluid channel is significantly improved compared to the smooth case, with the objective function value decreasing by approximately $37\%$, which aligns well with the simulation results.

According to the TO results in Fig. 3a, it is evident that the optimal topological structure varies with the inlet velocity. Here, we present more analysis in two cases: $\omega=0.1$ (mass transfer dominated; Fig. 5a) and $\omega=10$ (energy consumption dominated; Fig. 5b).

When $\omega=0.1$, the primary objective of the optimization process is to minimize $\sigma_{C}^{2}$, and we show the optimized concentration distributions for three different inlet velocities in Fig. 5a. We find that the vertical location of the solid baffle increases with $v$, which is consistent with the position where the concentration drops to 0 in Figs. 1b and d. This indicates that adding solid baffles enhances both the turbulent diffusion coefficient and the convective diffusion rate, which promotes the diffusion of solute toward regions of lower concentration, resulting in higher mass transfer rates and a more uniform concentration distribution. Another solid baffle is located at the lower wall, where the concentration boundary layer is more pronounced. As the thickness of the concentration boundary layer at this location decreases with the increase of $v$, the volume of the solid in this region also decreases accordingly.

For $\omega=10$, the primary objective of the TO is to reduce the $\Delta P$. The optimal topological structures (Fig. 5b) indicate that when the solid baffle is located at the bottom of the design domain $\Omega$, the pressure distribution of the process is similar to that of the smooth channel in Figs. 1b and d. This can be explained as that a relatively thick pressure boundary layer forms near the bottom wall due to the presence of inlet 2. When the solid is positioned near the bottom of the domain, its impact on the system’s pressure drop is relatively small. Since the position of the pressure boundary layer remains unchanged with variations in $v$, when the optimization process is dominated by minimizing $\Delta P$, the optimal topological structure is independent of $v$.

We compare the objective values before and after the TO in Fig. 5c. When the value of $\omega$ is relatively small, the process performance, controlled by reducing the concentration distribution variance, is significantly improved by the TO compared to the smooth channel. Moreover, for smaller $v$, the increase in system efficiency by the TO is more pronounced. The results for a larger value of the $\omega$ show that the external energy consumption of the smooth channel is lower than that of the optimized channel for inlet velocities $v<0.8\,\mathrm{m/s}$. However, for higher inlet velocities, the pressure distribution by the TO outperforms that of the smooth channel.

The turbulent mass transfer process (Fig. 1a) can be described by the Reynolds-averaged Navier-Stokes (RANS) equations [48, 49] with $i$ denoting $x$ and $y$ coordinates:

$$\frac{\partial\rho u_{i}}{\partial\mathbf{x}_{i}}=0,$$

(5)

$$u_{i}\frac{\partial\rho u_{j}}{\partial\mathbf{x}_{i}}=-\frac{\partial p}{\partial\mathbf{x}_{j}}+\frac{\partial}{\partial\mathbf{x}_{i}}\left(\left(\mu_{m}+\mu_{t}\right)\frac{\partial u_{j}}{\partial\mathbf{x}_{i}}\right)+f_{i},$$

(6)

$$u_{i}\frac{\partial\rho C}{\partial\mathbf{x}_{i}}=-\frac{\partial}{\partial\mathbf{x}_{i}}\left(\rho{(\gamma\cdot D}_{m}+D_{t})\frac{\partial C}{\partial\mathbf{x}_{i}}\right).$$

(7)

Eqs. (5)–(7) represent the continuity equation, momentum conservation equation, and species conservation equation, respectively. Here, $p$, $u$, and $\mathbf{x}$ denote pressure, velocity, and coordinate position, respectively. $C$ represents the time-averaged concentration. $\rho$, $\mu_{m}$, and $D_{m}$ are physical constants representing density, laminar viscosity, and laminar mass diffusivity, respectively. $\mu_{t}$ and $D_{t}$ represent turbulent viscosity and turbulent mass diffusivity, which can be obtained from the $k$-$\varepsilon$ two-equation model [50] and $\overline{c^{\prime 2}}$-$\varepsilon_{c^{\prime}}$ two-equation model [51], respectively. The two-equation models and their parameters can be found in Ref. [48]. Eqs. (5)–(6) and the closure equation for turbulent viscosity $\mu_{t}$ constitute the CFD equations used to calculate the velocity and pressure distributions. The CMT equations include the CFD equations, Eq. (7) for concentration calculation, along with the closure equations for turbulent mass diffusivity $D_{t}$.

Additionally, to account for the influence of changing the system’s topological structure on the physical fields, we introduce a source term $f_{i}$ in the momentum conservation equation (Eq. (6)) to represent the resistance effect of solid baffles on the flow field. Moreover, in the component conservation equation, a fluid volume fraction $\gamma$ within the range [0, 1] is introduced to represent the effect of the solid baffle on the component diffusion coefficient. Specifically, $f_{i}$ represents the frictional drag generated by solid baffles and can be solved using a density model [19]:

$$f_{i}=-\alpha u_{i},$$

where $\alpha$ is the interpolation function of the porosity $\gamma$:

$$\alpha=\alpha_{\text{min}}+(\alpha_{\text{max}}-\alpha_{\text{min}})q\frac{1-\gamma}{1+\gamma}.$$

$\alpha$ represents the solid density, and $q$ is a positive integer used to adjust the shape of the interpolation curve. In this study, $q$ and $\alpha_{\text{min}}$ are set to 1 and 0, respectively [52]. $\alpha_{\text{max}}$ is a large value that approximates the fluid-to-solid transition, set to 600,000. When $\gamma$ is 0, the corresponding fluid density is $\alpha_{\text{max}}$, resulting in significant frictional resistance $f_{i}$, making the domain effectively solid-like. Conversely, when $\gamma$ is 1, the corresponding region is fluid. In this study, we set the velocity inlet and pressure outlet (with static pressure fixed at 0) [24]. More computational details can be found in Ref. [24].

To train the neural operators, various topological structures $\gamma(x,y)$ are randomly generated and evenly assigned to 9 different inlet velocities ($v\in\{0.1,0.2,\dots,0.9\}\,\text{{m/s}}$). To obtain the corresponding concentration $C$ and pressure $p$ distributions for the given $v$ and $\gamma(x,y)$, the mechanism model (Section 4.1) is solved in the commercial CFD software package FLUENT $14.5^{TM}$. The value of $\gamma(x,y)$ is stored using user-defined memory, while $\mu_{t}$, $D_{t}$, and $f_{i}$ in Eqs. (6)–(7) are implemented through user-defined functions. The total pressure $p$, including static and dynamic pressures, is normalized as $P=\frac{p}{2000}+0.26$. Each data point in the dataset includes the coordinates $(x,y)$, along with their corresponding values of $(\gamma,P,C)$. The initial dataset is categorized into training, testing, and interpolation/extrapolation datasets based on the value of the inlet velocity $v$. More details on the initial dataset are provided in Section S1.

Fourier-DeepONet, illustrated in Fig. 2 Step 1, was developed based on three neural network frameworks: the deep operator network (DeepONet) [31], the Fourier neural operator (FNO) [53], and U-FNO (a block combining Fourier neural operator and U-Net) [54]. DeepONet is used to encode the input variables and functions, while FNO and U-FNO apply Fourier transforms to decode DeepONet’s output to obtain the model output. Fourier-DeepONet exhibits superior generalizability and better accuracy in learning neural operators in high-dimensional spaces [45].

The vanilla DeepONet consists of two components: a branch net and a trunk net used for encoding the inputs. In this study, the trunk net takes the inlet velocity $v$ as input, while the branch net takes the distribution of fluid volume $\gamma(x,y)$ of the design domain $\Omega$ (Fig. 1a) as input. Since the solutions of the mechanism model are used as data for training Fourier-DeepONet, the dimension of the input function $\gamma(x,y)$ corresponds to the number of grid nodes of the mechanism model, which is 121 (mesh nodes in the $x$-direction) $\times$ 101 (mesh nodes in the $y$-direction). The outputs of the branch net and the trunk net are denoted as $\mathbf{b}$ and $\mathbf{t}$, respectively:

$$\mathbf{b}=B[P\left\langle\gamma(x,y)\right\rangle]\in\mathbb{R}^{L_{0}\times H_{0}\times C},$$

$$\mathbf{t}=T(v)\in\mathbb{R}^{C},$$

where $B$ and $T$ represent two linear transformations, and $P\left\langle\cdot\right\rangle$ denotes the padding operation. $L_{0}\times H_{0}$ represents the output dimensionality of a single channel in the branch net, where in this study $L_{0}=128$ and $H_{0}=112$. $C$ represents the number of channels, which corresponds to the width of the operator layers. Additionally, as shown in Fig. 2, DeepONet computes $\mathbf{z}_{0}$ by combining the outputs $\mathbf{b}$ and $\mathbf{t}$ through element-wise multiplication:

$$\mathbf{z}_{0}=\mathbf{b}\odot\mathbf{t}.$$

We then utilize an FNO layer and a U-FNO layer to decode the output of DeepONet. The outputs of the FNO layer, $\mathbf{z}_{1}$, and the two U-FNO layers, $\mathbf{z}_{2}$ and $\mathbf{z}_{3}$, are computed as

$$\mathbf{z}_{1}=\sigma(\mathcal{F}^{-1}\left(R_{1}\cdot\mathcal{F}\left(\mathbf{z}_{0}\right)\right)+W_{1}\mathbf{z}_{0}+b_{1}),$$

$$\mathbf{z}_{2}=\sigma(W_{2}^{\prime}(\mathcal{F}^{-1}\left(R_{2}\cdot\mathcal{F}\left(\mathbf{z}_{1}\right)\right)+\mathcal{U}_{2}\left(\mathbf{z}_{1}\right)+W_{2}\mathbf{z}_{1}+b_{2})),$$

$$\mathbf{z}_{3}=\sigma(W_{3}^{\prime}(\mathcal{F}^{-1}\left(R_{3}\cdot\mathcal{F}\left(\mathbf{z}_{2}\right)\right)+\mathcal{U}_{3}\left(\mathbf{z}_{2}\right)+W_{3}\mathbf{z}_{2}+b_{3})).$$

In the FNO and U-FNO layers, $\mathcal{F}$ represents the two-dimensional Fast Fourier Transform (FFT), $\mathcal{F}^{-1}$ denotes the inverse two-dimensional FFT, $\sigma$ is the activation function, $W_{i}$ are weight matrices, $R_{i}$ are complex valued tensors, and $b_{i}$ are bias. In the U-FNO layer, $\mathcal{U}_{i}$ denotes a U-Net layer. Since the input and output dimensions of the Fourier and U-Net layers must be consistent, a linear layer is added before the activation function in the U-FNO layer to transform the output function to the same dimension as the concentration and pressure distributions. $W^{\prime}_{i}$ represents the weight matrix of this linear layer.

The projection layer in the decoding part is used to perform nonlinear transformations and slicing operations on the decoded output function $\mathbf{z}_{3}$ from the Fourier layer, resulting in the model’s output, i.e., the physical fields $U(x,y)=C$ or $P$. The projection layer can be expressed as

$$U\left(x,y\right)=S\left\langle(W_{P2}\sigma(W_{P1}\mathbf{z}_{3}+b_{P1})+b_{P2})\right\rangle,$$

where $S\left\langle\cdot\right\rangle$ denotes the slicing operation, $W_{P1}$ and $W_{P2}$ are the weight matrices, and $b_{P1}$ and $b_{P2}$ are biases of the projection layer. The dimensions of the concentration distribution $C(x,y)$ and the pressure distribution $P(x,y)$ are consistent with the number of nodes in the CFD simulation, i.e., 201 (number of nodes in the $x$-direction) $\times$ 101 (number of nodes in the $y$-direction).

As shown in Fig. 1c, the pressure distribution is more sensitive to changes in the topological structure than the concentration distribution, which requires more channels for $\mathcal{G}_{P}$ than for $\mathcal{G}_{C}$. In this study, the number of channels used for $\mathcal{G}_{C}$ and $\mathcal{G}_{P}$ was set to 32 and 48, respectively. For other model hyperparameters of DeepONet and the projection layer $Q$, please refer to Refs. [31, 45], and for other model hyperparameters of FNO and U-FNO, refer to Refs. [54, 53]. The training loss trajectories of $\mathcal{G}_{C}$ and $\mathcal{G}_{P}$ are shown in Section S3 Fig. S2.

As shown in Fig. 2 Step 2, a topological structure is constructed by a neural network, which consists of two layers, each containing 140 neurons. $\gamma_{0}$ represents the preliminary predicted porosity values by the neural network. As the porosity $\gamma$ must take values of 0 or 1 for any input $(v,x,y)$, to obtain a reasonable topological structure, the preliminary output $\gamma_{0}$ needs to be transformed to 0/1. In this study, the sigmoid function is used for the 0/1 transformation, and the topological structure is computed as

$$\gamma(x,y)=sigmoid(\alpha\cdot(\gamma_{0}(x,y)-\bar{\gamma_{0}})),$$

where $\bar{\gamma_{0}}$ represents the average value of $\gamma_{0}$ at different coordinate points $(x,y)$ for the same velocity boundary. $\alpha$ is a hyperparameter, and a larger value of $\alpha$ indicates a sharper change between the output of 0 and 1. In this study, $\alpha=10$.

The trained $\mathcal{G}_{P}$ and $\mathcal{G}_{C}$ neural operators take the $v$ and the topological structure as the inputs to predict the pressure and concentration distributions, respectively. Then $\sigma_{C}^{2}$, $\Delta P$, and $Obj$ are computed according to Eqs. (2)–(4). The training losses $\mathcal{L}_{\sigma_{C}^{2}}$, $\mathcal{L}_{\Delta P}$, and $\mathcal{L}_{obj}$ are computed from $\sigma_{C}^{2}$, $\Delta P$, and $Obj$ under different inlet velocities:

$$\mathcal{L}_{\sigma_{C}^{2}}=\frac{1}{N}\sum_{n=1}^{N}{{\sigma_{C}^{2}}^{(n)}},$$

$$\mathcal{L}_{\Delta P}=\frac{1}{N}\sum_{n=1}^{N}{{\Delta P}^{(n)}},$$

$$\mathcal{L}_{obj}=\frac{1}{N}\sum_{n=1}^{N}{Obj^{(n)}},$$

where the superscript $n$ represents different data points in the training dataset. The total inequality constraint of the solid volume $\mathcal{L}_{vol}$ is computed based on Eq. (1):

$$\mathcal{L}_{vol}=\frac{1}{N}\sum_{n=1}^{N}\min\left\{\iint_{\Omega}\left(1-\gamma^{(n)}\right)dV-0.1\times V_{\Omega},0\right\}.$$

The optimization of a neural topology does not require observational data, and the training loss $\mathcal{L}$ is composed of topological structure constraints and the objective function:

$$\mathcal{L}=\mathcal{L}_{obj}+\lambda\cdot\mathcal{L}_{vol},$$

where $\lambda$ is a weight coefficient and is set to $1\times 10^{3}$ to ensure the optimal topology satisfies the volume constraint.

During the training, the initial learning rate decays from $1\times 10^{-3}$ to $1\times 10^{-4}$ in 5000 steps. The activation function and optimizer are ReLU and Adam, respectively. When $\omega=1$, the training error curve of the TO framework is shown in Fig. S3. We observe that $\mathcal{L}_{vol}$ approaches 0, indicating that the TO results satisfy the inequality constraints of solid volume.

After performing Step 1 and Step 2 using the initial dataset in Section S1, the TO result I is obtained, and some optimized structures are depicted in Fig. S5a. In Step 2, we use neural operators to replace the mechanism model for enabling fast prediction of the objectives. Hence, a good prediction accuracy of neural operators is crucial, and the errors are listed in Table S1. $\mathcal{G}_{P}$ shows lower prediction accuracy for the optimized topological structure compared to $\mathcal{G}_{C}$, because the training dataset size of $\mathcal{G}_{P}$ is equal to the number of topologies in the dataset, regardless of the inlet velocities, and thus it is much smaller than the dataset size of $\mathcal{G}_{C}$. Moreover, as suggested by the TO result I, the optimal structure would only have one or two solid structures near the bottom or left boundaries of the design domain. However, in our randomly-generated initial dataset, the structures do not satisfy this pattern.

To improve the prediction accuracy of $\mathcal{G}_{P}$, especially for structures similar to the optimal structure, more data are generated as follows. In the first round of data augmentation, the solid positions are chosen the same as those in Fig. S5a, while the solid shapes are randomly generated. We generate ten new structures for $v=0.1,0.2,\dots,0.9\,\mathrm{m/s}$, and we show three examples in Fig. S4. To reduce the computational cost, we only use these structures to generate the training data of $\mathcal{G}_{P}$. Then, TO result II is obtained through transfer learning of the $\mathcal{G}_{P}$ neural operator and the neural topology.

In the second round of data augmentation, we use the 27 optimal structures obtained from TO result II to generate new training data for both $\mathcal{G}_{C}$ and $\mathcal{G}_{P}$.