|
| 1 | +--- |
| 2 | +title: Adjoint design optimization of a turbulent 3D pipe bend |
| 3 | +permalink: /tutorials/Inc_Turbulent_Bend_Opt/ |
| 4 | +written_by: Nijso Beishuizen |
| 5 | +for_version: 8.0.1 |
| 6 | +revised_by: |
| 7 | +revision_date: |
| 8 | +revised_version: |
| 9 | +solver: INC_RANS |
| 10 | +requires: SU2_CFD, FADO |
| 11 | +complexity: advanced |
| 12 | +follows: Inc_Turbulent_Bend |
| 13 | +--- |
| 14 | + |
| 15 | + |
| 16 | + |
| 17 | + |
| 18 | +#### Figure (1): The 90 degree bend with the FFD box. |
| 19 | + |
| 20 | + |
| 21 | +## Goals |
| 22 | + |
| 23 | +This tutorial closely follows the design optimization setup of the 2D mixing channel, see the tutorial on design optimization for [Species Transport](/tutorials/Species_Transport/). We will use the python framework [FADO](https://github.com/su2code/FADO) to set up the optimization problem. In this tutorial however we will optimize the pressure drop of the 90 degree pipe bend. The CFD results were already discussed previously in [this tutorial](/tutorials/Inc_Turbulent_Bend/). Here we focus on the following aspects: |
| 24 | +* Setup of the FFD box for a 3D problem using SU2_DEF. |
| 25 | +* Setup of the FADO problem: |
| 26 | + - using previous solutions of the CFD and adjoint solver to reduce computing time. |
| 27 | + - automatically add the correct keywords for the FFD box (type and degrees of freedom). |
| 28 | + - discussion of mesh quality issues and improvement using FFD constraints. |
| 29 | +* Presentation of the final setup and results. |
| 30 | + |
| 31 | +If you have not done so already, please first have a look at the prerequisite tutorials for a better understanding. |
| 32 | + |
| 33 | +## Prerequisites |
| 34 | + |
| 35 | +Besides the python library FADO, you will need to compile SU2 with automatic differentiation support. Note that the script provided uses 8 cores, so you need to compile with mpi support as well as enable autodiff: |
| 36 | +``` |
| 37 | +./meson.py build --optimization=2 -Ddebug=false -Denable-autodiff=true -Dwith-mpi=enabled --prefix=/home/user/Codes/su2_github_develop/su2/ |
| 38 | +``` |
| 39 | + |
| 40 | +## Resources |
| 41 | + |
| 42 | +You can find the resources for this tutorial in the folder [design/Inc_Turbulent_Bend_Wallfunctions](https://github.com/su2code/Tutorials/tree/master/design/Inc_Turbulent_Bend_Wallfunctions) and the respective subfolders. Note that the setup used in this tutorial uses a pipe bend with a shorter pipe length before and after the bend to reduce the computing time. We have also merged all wall boundaries and all symmetry planes into one wall boundary and one symmetry plane. The gmsh file is provided in the repository so you can create your own pipe bend(s) with it. |
| 43 | + |
| 44 | +## 1. Basic setup of FFD box and FADO |
| 45 | +Usually, designs are created with a CAD tool. These designs are then discretized into a computational mesh for CFD analysis. When optimizing a design, we then only have the discretized mesh available. We could manipulate the mesh nodes directly but this does not lead to very smooth deformations. Instead we modify our mesh using FFD boxes. The nodes of the FFD box are moved according to the design sensitivities, and the mesh nodes inside the FFD box are then smoothly deformed using Bezier curves (default) or B-splines. |
| 46 | +Figure (1) above shows the setup that we will be using. |
| 47 | + |
| 48 | +### Creation of the FFD box |
| 49 | +The FFD box can be created and added to the .su2 mesh file using SU2_DEF. The important parameters to add to the configuration file are: |
| 50 | + |
| 51 | +``` |
| 52 | +FFD_DEFINITION= (BOX, \ |
| 53 | +-0.06, -0.009, -0.001, \ |
| 54 | + 0.209, -0.009, -0.001, \ |
| 55 | + 0.209, 0.060, -0.001, \ |
| 56 | +-0.06, 0.060, -0.001, \ |
| 57 | +-0.06, -0.009 ,0.27, \ |
| 58 | + 0.209 ,-0.009, 0.27, \ |
| 59 | + 0.209, 0.060, 0.27, \ |
| 60 | +-0.06, 0.060, 0.27 ) |
| 61 | +
|
| 62 | +FFD_DEGREE= (5,5,5) |
| 63 | +DV_KIND= FFD_SETTING |
| 64 | +DV_MARKER= ( wall, symmetry) |
| 65 | +DV_PARAM= (1.0) |
| 66 | +``` |
| 67 | +Our FFD box is 5x5x5 cells, or 6x6x6=216 nodes. With each 3 degrees of freedom for the x-,y- and z-direction, we get a total of 648 d.o.f. The box is slightly larger than our original bend, but most importantly the symmetry plane is completely inside the FFD box. |
| 68 | +We run the command |
| 69 | + |
| 70 | +``` |
| 71 | +$ SU2_DEF sudo_0_add_FFD_box.cfg |
| 72 | +``` |
| 73 | + |
| 74 | +We will only use the file *sudo_0_add_FFD_box.cfg* to create the FFD box. This command will create a new .su2 mesh called *mesh_out.su2* that has the definition of the FFD box added to the file. Note that at this stage we need to provide the boundaries inside the FFD box that are allowed to deform using the keyword **DV_MARKER**. The nodes on the boundaries inside the FFD box are part of the information that is added with the FFD box to the mesh_out.su2 mesh file. |
| 75 | + |
| 76 | + |
| 77 | +### Setup the FADO script optimization.py |
| 78 | + |
| 79 | +Previously the python script *set_ffd_design_var.py* was used in the tutorial on the optimization of the mixing channel. We can however create the correct entries with a couple of simple python commands and add them directly to the main python script that FADO will use, *optimization.py*. We already used *optimization.py* before in the tutorial [Species Transport](/tutorials/Species_Transport/) and we will highlight some changes here. |
| 80 | + |
| 81 | +For the optimization, we need to modify 2 things in our config file: **DV_KIND** and **DV_PARAM**. The keyword **DV_PARAM** contains N entries of the form *( BOX, 0, 0, 0, 1.0, 0.0, 0.0 );* |
| 82 | + Note that the meaning of the entries are *( FFD_BoxTag, i_Ind, j_Ind, k_Ind, x_Disp, y_Disp, z_Disp )* , meaning that after the keyword **BOX**, we get the 3 indices i,j,k of the FFD box, followed by the allowed displacement of that index in the x-,y- and z-direction. Mesh nodes in the symmetry plane only move in the symmetry plane, so symmetry is preserved. |
| 83 | + |
| 84 | +The list of FFD control points can be created using: |
| 85 | +```python |
| 86 | +s = "FFD_CONTROL_POINT" |
| 87 | +ffd_string = s |
| 88 | +for i in range((DX**NDIM)*NDIM - 1): |
| 89 | + ffd_string = ffd_string + ", " + s |
| 90 | +``` |
| 91 | +In the sudo.cfg file we have a placeholder which has the form: |
| 92 | +``` |
| 93 | +DV_KIND= __FFD_CTRL_PTS__ |
| 94 | +``` |
| 95 | +And we automatically replace the placeholder **\_\_FFD_CTRL_PTS\_\_** in the fado script during run-time using the command |
| 96 | +``` |
| 97 | +replace_dv_kind = Parameter([ffd_string], LabelReplacer("__FFD_CTRL_PTS__")) |
| 98 | +``` |
| 99 | + |
| 100 | +For the **DV_PARAM** string, we can use a similar piece of python code: |
| 101 | + |
| 102 | +```python |
| 103 | +dv_param_string="" |
| 104 | +for idim in range(NDIM): |
| 105 | + xdim = ydim = zdim = "0.0" |
| 106 | + if (idim==0): xdim="1.0" |
| 107 | + elif (idim==1): ydim="1.0" |
| 108 | + elif (idim==2): zdim="1.0" |
| 109 | + for k in range(DX): |
| 110 | + for j in range(DX): |
| 111 | + for i in range(DX): |
| 112 | + s = "( BOX, " + str(i) + ", " + str(j) + ", " + str(k) + ", " + xdim + ", " + ydim + ", " + zdim + " );" |
| 113 | + dv_param_string += s |
| 114 | +``` |
| 115 | +which is replaced using: |
| 116 | +```python |
| 117 | +replace_dv_param =Parameter([dv_param_string], LabelReplacer("__FFD_PARAM__")) |
| 118 | +``` |
| 119 | + |
| 120 | +The only thing we need to take care of now is to define the correct number of FFD nodes DX=6 and the correct dimension of the problem, NDIM=3. You could even get these values from the .su2 mesh file if you want! |
| 121 | + |
| 122 | +In our config file sudo.cfg, we also use the following setting: |
| 123 | +``` |
| 124 | +FFD_CONTINUITY= 1ST_DERIVATIVE |
| 125 | +``` |
| 126 | +Which means that when the sides of the FFD box cuts through the mesh, then at the cut, we want the mesh deformation to stay first order continuous. If you do not do this, you might get sharp jumps in your mesh at the location of the intersection. |
| 127 | + |
| 128 | +Unfortunately, if we use this unconstrained setup, the mesh deformation will be as shown in Figure (2) below: |
| 129 | + |
| 130 | + |
| 131 | +#### Figure (2): The unconstrained deformation leads to mesh cells that have collapsed on the symmetry plane. This leads to convergence issues at the next design iteration. |
| 132 | + |
| 133 | +As you can see, the row of cells just above the symmetry plane has collapsed onto the symmetry plane, leading to a very bad mesh. In the next design iteration, the simulations do not converge anymore and the optimization procedure stops. So we need some additional constraints on the movement of the nodes of the FFD box. |
| 134 | + |
| 135 | +### constrained FFD deformation |
| 136 | + |
| 137 | +The problem here is that the mesh nodes on the symmetry plane are not allowed to move vertically, they only move horizontally outward inside the symmetry plane. Unfortunately, the FFD deformation is such that mesh nodes just above the symmetry plane are moved down, almost on the symmetry plane. To improve the quality of the mesh deformation, we disallow the vertical movement of the FFD box nodes on the nodes in the bottom plane of the FFD box, with j-index 0 and vertical displacement (0,1,0). In Figure (3), the plane with index j=0 is the bottom plane, indicated in yellow. So we remove entries in the **DV_PARAM** list of the form *(BOX, i_Ind, 0, k_Ind, 0,1,0)*. Additionally, we also disallow the vertical movement of the nodes in the plane j=1. |
| 138 | +Since we disallow vertical movement in 2x(6x6)=72 nodes in the planes with $$j=0$ and $j=1$$, The total degrees of freedom is then $$648 - 72 = 576$$ d.o.f. |
| 139 | + |
| 140 | + |
| 141 | + |
| 142 | +#### Figure (3): The 90 degree bend with the FFD box. |
| 143 | + |
| 144 | + |
| 145 | +The number of design variables needs to be reduced by 2x6x6 and we need to remove 2x6x6 strings from the ffd_string: |
| 146 | +```python |
| 147 | +nDV = nDV - 2*DX*DX |
| 148 | +ffd_string.replace(s+", ","",2*DX*DX) |
| 149 | +``` |
| 150 | + |
| 151 | +Additionally, we need to remove the vertical movement of the nodes in the plane with j=[0,1]: |
| 152 | + |
| 153 | +```python |
| 154 | +jlist = [0,1] |
| 155 | +dof = "0.0, 1.0, 0.0" |
| 156 | + |
| 157 | +for j in jlist: |
| 158 | + for k in range(DX): |
| 159 | + for i in range(DX): |
| 160 | + remove_dof = "( BOX, " + str(i) + ", " + str(j) + ", " + str(k) + ", " + dof + " )" |
| 161 | + print("removing ", remove_dof) |
| 162 | + dv_param_string = dv_param_string.replace(remove_dof+";", "", 1) |
| 163 | + # in case the plane is at the end, the string does not have a semicolon |
| 164 | + dv_param_string = dv_param_string.replace(remove_dof, "", 1) |
| 165 | +``` |
| 166 | + |
| 167 | +Another modification to the configuration file that we would like to add is to restart from the solution of the previous design. SU2 will automatically interpolate the solution to the nearest neighbor if the mesh coordinates do not match. Since we might not have an initial solution to start with, we would like to switch from **RESTART= NO** to **RESTART= YES** after the first iteration. |
| 168 | + |
| 169 | +We can do this with a simple *sed* command that searches and replaces the string in config.cfg in the working directory and then copies the config file back to the base dir: |
| 170 | + |
| 171 | +```python |
| 172 | +restart_yes="sed -i 's/RESTART_SOL= NO/RESTART_SOL= YES/' config.cfg && cp " + configCopy + " ../../" |
| 173 | +``` |
| 174 | + |
| 175 | +The easiest way to perform the update is to simply add it as another ExternalRun process: |
| 176 | +```python |
| 177 | +update_restart = ExternalRun("UPDATE_RESTART",restart_yes,False) # True means sym links are used for addData |
| 178 | +update_restart.addData(configCopy) |
| 179 | +``` |
| 180 | +This means that we will call this function every design iteration, but only in the first design iteration will it have an effect. |
| 181 | + |
| 182 | +FADO copies the restart file and all other necessary files from the base directory (where you start your script) into the working directory. Here, the working directory is called *OPTIM*. Inside *OPTIM* we have subdirectories for the different runs, i.e. *OPTIM/DEFORM*, OPTIM/DIRECT*, *OPTIM/ADJOINT* and *OPTIM/DOT*. When all runs in this directory are done, they are archived in *DSN_001*, *DSN_002*, etc. |
| 183 | + |
| 184 | +So what we will do now is every time after we run the primal solver, we copy the restart file from the working directory *OPTIM/DIRECT* back to the base dir. |
| 185 | + |
| 186 | +```python |
| 187 | +cfd_command = "mpirun -n " + ncores + " " + su2_run + "SU2_CFD " + configMaster + "&& cp restart.csv ../../solution.csv" |
| 188 | +``` |
| 189 | + |
| 190 | +Note that the primal solver by default saves the file restart.csv (keyword: **RESTART_FILENAME= restart**), but reads the file solution.csv (keyword: **SOLUTION_FILENAME= solution**). Note that we save ASCII files here, which is determined by the option **OUTPUT_FILES= RESTART_ASCII**. When reading ASCII restart files, we also need **READ_BINARY_RESTART= NO**. |
| 191 | +Please go through the *optimization.py* script now and check the settings. By default, the above procedure to restart from a previous solution is not active. You can activate it now if you wish. |
| 192 | + |
| 193 | +## 2. FADO optimized result |
| 194 | + |
| 195 | +We run the optimization script: |
| 196 | + |
| 197 | +``` |
| 198 | +$ python optimization.py |
| 199 | +``` |
| 200 | + |
| 201 | +We have set it to run 5 design iterations. Please note that the optimization was set up for a parallel run. |
| 202 | + |
| 203 | +The objective values for the 5 design iterations are given in the table below. |
| 204 | + |
| 205 | + |
| 206 | +| iteration | $$\Delta P$$, total| $$\Delta P$$, bend| gain, total| gain, bend | |
| 207 | +| --------|--------|--------|--------|--------| |
| 208 | +|0|89.0 [Pa]|20.7 [Pa]| -| - | |
| 209 | +|1|84.9 [Pa]|16.6 [Pa]|4.6 \%|19.8 \%| |
| 210 | +|2|82.1 [Pa]|13.8 [Pa]|7.8 \%|33.3 \%| |
| 211 | +|3|81.6 [Pa]|13.3 [Pa]|8.3 \%|35.7 \%| |
| 212 | +|4|79.9 [Pa]|11.6 [Pa]|10 \%|44 \%| |
| 213 | + |
| 214 | + |
| 215 | +We see that the global pressure drop between the inlet and the outlet reduces from 89 Pa to 80 Pa, a reduction of more than 10 \%. However, since we are optimizing only the bend, we should subtract the pressure drop of the straight parts and only consider the pressure drop of the bend for a more fair comparison. In paraview, we can integrate the pressure in a 2D slice at the start and end of the bend. The pressure drop of the bend comes to $$\Delta P = 20.7 Pa$$. That means that the reduction of pressure drop in the bend is actually 44 \% ! |
| 216 | + |
| 217 | + |
| 218 | +#### Figure (4): Optimized bend after 5 design iterations, with the deformed FFD box. |
0 commit comments