Skip to content

Commit 5f6fd9b

Browse files
rsanferrsanfer
committed
Added adjoint FSI with python.
1 parent 9a73722 commit 5f6fd9b

File tree

1 file changed

+378
-0
lines changed

1 file changed

+378
-0
lines changed

‎_tutorials/multiphysics/Adjoint_FSI_Python.md‎

Lines changed: 378 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,378 @@
1+
---
2+
title: Static Fluid-Structure Interaction (FSI) Adjoint via Python-wrapper
3+
permalink: /tutorials/Adjoint_FSI_Python/
4+
written_by: rsanfer
5+
for_version: 7.0.3
6+
revised_by: rsanfer
7+
revision_date: 2020-03-04
8+
revised_version: 7.0.3
9+
solver: Multiphysics
10+
requires: pysu2ad library
11+
complexity: advanced
12+
follows: Static_FSI
13+
follows2: Static_FSI_Python
14+
userguide: Build-SU2-Linux-MacOS
15+
---
16+
17+
This tutorial uses SU2's python wrapper and its native adjoint solvers for incompressible flow and solid mechanics to solve a steady-state, **adjoint** Fluid-Structure Interaction problem. This document will cover:
18+
- Operating with the ad version of the pysu2 library
19+
- Extracting the adjoints of the flow loads and structural displacements from two different python instances of SU2
20+
- Exchanging adjoint information between the two instances
21+
22+
In this tutorial, we will solve the adjoint of the problem presented in the [Static FSI with Python](../Static_FSI_Python/) tutorial, which is a pre-requisite in order to be able to run this tutorial.
23+
24+
### Resources
25+
26+
You can find the resources for this tutorial in the folder [python_fsi](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi) of the [Tutorials repository](https://github.com/su2code/Tutorials). There is a [python script](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/run_fsi_adjoint.py) and two sub-config files for the [flow AD subproblem](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/config_channel_ad.cfg) and [structural AD subproblem](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/config_cantilever_ad.cfg).
27+
28+
Moreover, you will need to reuse the two mesh files from the [Static FSI with Python](../Static_FSI_Python/) tutorial and the solution files generated when running it.
29+
30+
### Background
31+
32+
For this tutorial, you will need to use advanced features of SU2, in particular the python-wrapped version of the AD code, which needs to be built from source.
33+
34+
This tutorial has been tested on a linux system with the following specs
35+
36+
```
37+
Linux kernel 5.3.18-1-MANJARO
38+
GCC compilers version 9.2.0
39+
Open MPI version 4.0.2
40+
Python version 3.8.1
41+
SWIG version 4.0.1
42+
CoDiPack version 1.8
43+
```
44+
45+
compiling the code from source using the following meson settings
46+
47+
```
48+
./meson.py build -Dwith-mpi=enabled -Denable-autodiff=true -Denable-pywrapper=true
49+
```
50+
51+
It requires an adequate setup of the system and a correct linkage of both the ```mpi4py``` and ```swig``` libraries. For questions, updates, or notes on the usability of this tutorials on different systems or configurations, please use the comment section below.
52+
53+
#### Mesh Description
54+
55+
The cantilever is discretized using 1000 4-node quad elements with linear interpolation. The fluid domain is discretized using a finite volume mesh with 7912 nodes and 15292 triangular volumes. The wet interface is matching between the fluid and structural domains.
56+
57+
#### Configuration File Options
58+
59+
We reuse the [flow](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/config_channel_ad.cfg) and [structural](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/config_cantilever_ad.cfg) config files from the [Static FSI with Python](../Static_FSI_Python/) tutorial. However, it is necessary to add some changes to run the discrete adjoint solver.
60+
61+
First, we need to enable adjoint mode on the flow config file
62+
63+
```
64+
MATH_PROBLEM = DISCRETE_ADJOINT
65+
OBJECTIVE_FUNCTION = DRAG
66+
```
67+
68+
and the structural config file
69+
70+
```
71+
MATH_PROBLEM = DISCRETE_ADJOINT
72+
```
73+
74+
where the drag coefficient is defined as the functional for which we want to compute the adjoint.
75+
76+
Next, and same as for the primal tutorial, SU2 will see each instance as a single-zone problem. Therefore, it is necessary to set the output appropriately in order to prevent overwriting files. For the flow AD config
77+
78+
```
79+
OUTPUT_FILES = (RESTART, PARAVIEW)
80+
SOLUTION_FILENAME = solution_fsi_steady_0
81+
RESTART_FILENAME = restart_fsi_steady_0
82+
83+
SOLUTION_ADJ_FILENAME = solution_ad_fsi_steady_0
84+
RESTART_ADJ_FILENAME = restart_ad_fsi_steady_0
85+
VOLUME_ADJ_FILENAME = fsi_ad_steady_0
86+
```
87+
88+
and for the structural config
89+
90+
```
91+
OUTPUT_FILES = (RESTART, PARAVIEW)
92+
SOLUTION_FILENAME = solution_fsi_steady_1
93+
RESTART_FILENAME = restart_fsi_steady_1
94+
95+
SOLUTION_ADJ_FILENAME = solution_ad_fsi_steady_1
96+
RESTART_ADJ_FILENAME = restart_ad_fsi_steady_1
97+
VOLUME_ADJ_FILENAME = fsi_ad_steady_1
98+
```
99+
100+
Next, we choose the output for the history files. In this tutorial, we are interested on the derivative of the drag coefficient with respect to the Young's modulus of the cantilever, and therefore for the flow simulation it is enough to retrieve the residuals on the history file
101+
102+
```
103+
HISTORY_OUTPUT = ITER, RMS_RES
104+
CONV_FILENAME= history_ad_0
105+
```
106+
107+
while for the structural config we need to request the sensitivities as well
108+
109+
```
110+
HISTORY_OUTPUT = ITER, RMS_RES, SENSITIVITY
111+
CONV_FILENAME= history_ad_1
112+
```
113+
114+
Also, we monitor in this case adjoint quantities for the flow simulation
115+
116+
```
117+
CONV_FIELD = RMS_ADJ_PRESSURE, RMS_ADJ_VELOCITY-X, RMS_ADJ_VELOCITY-Y
118+
CONV_RESIDUAL_MINVAL = -12
119+
```
120+
121+
and the structural simulation
122+
123+
```
124+
CONV_FIELD = ADJOINT_DISP_X, ADJOINT_DISP_Y
125+
CONV_STARTITER = 2
126+
CONV_RESIDUAL_MINVAL = -7
127+
```
128+
129+
Finally, and the same as for the primal case, it is necessary to indicate to the structural solver that the boundary solution will be used to update a fluid field. This is done by setting the config options
130+
131+
```
132+
MARKER_FLUID_LOAD = ( feabound )
133+
MARKER_DEFORM_MESH = ( feabound )
134+
```
135+
136+
on the structural config file.
137+
138+
#### Applying coupling conditions to the individual domains
139+
140+
As usual for adjoint problems, the first step is to rename the restart files from the primal run as solution files, ```restart_fsi_steady_0.dat``` → ```solution_fsi_steady_0.dat``` and ```restart_fsi_steady_1.dat``` → ```solution_fsi_steady_1.dat```.
141+
142+
The key part of this tutorial is the [python script](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/run_fsi_adjoint.py) to run the adjoint FSI problem. Please take a moment to evaluate its contents, as we will go through some of its most important aspects.
143+
144+
First, we will need to import the SU2 adjoint library along with mpi4py, using
145+
146+
```
147+
import pysu2ad as pysu2
148+
from mpi4py import MPI
149+
```
150+
151+
We define the names of the config files required by SU2 using
152+
153+
```
154+
flow_filename = "config_channel_ad.cfg"
155+
fea_filename = "config_cantilever_ad.cfg"
156+
```
157+
158+
We will exemplify the initialization of SU2 using the flow domain. First, we create a single-zone adjoint driver object using
159+
160+
```
161+
FlowDriver = pysu2.CDiscAdjSinglezoneDriver(flow_filename, 1, comm);
162+
```
163+
164+
which is analogous to the SU2 driver in the C++ executable. We use a ```comm``` that is imported from ```mpi4py```
165+
166+
```
167+
comm = MPI.COMM_WORLD
168+
```
169+
170+
and the ```flow_filename``` variable previously defined. Next, identifyin the FSI boundary is analogous to the primal script and will not be repeated here.
171+
172+
Now, the major differences with respect to the primal case are presented. First, we need to initialize the cross dependency that is applied as a source term into the flow domain to 0, using
173+
174+
```
175+
fea_sensitivities=[]
176+
for iVertex in range(nVertex_Marker_Flow):
177+
fea_sensitivities.append([0.0, 0.0, 0.0])
178+
```
179+
180+
We start the FSI loop and limit it to 15 iterations
181+
182+
```
183+
for i in range(15):
184+
```
185+
186+
and the source term corresponding to the flow load adjoint is applied to the fluid domain. In the first iteration, this will be a zero-vector, but that will not be the case for subsequent iterations.
187+
188+
```
189+
FlowDriver.SetFlowLoad_Adjoint(FlowMarkerID, 0, fea_sensitivities[1][0], fea_sensitivities[1][1], fea_sensitivities[1][2])
190+
FlowDriver.SetFlowLoad_Adjoint(FlowMarkerID, 1, fea_sensitivities[0][0], fea_sensitivities[0][1], fea_sensitivities[0][2])
191+
for iVertex in range(2, nVertex_Marker_Flow):
192+
FlowDriver.SetFlowLoad_Adjoint(FlowMarkerID,iVertex,fea_sensitivities[iVertex][0],fea_sensitivities[iVertex][1],fea_sensitivities[iVertex][2])
193+
```
194+
195+
The flow adjoint iteration is run now using
196+
197+
```
198+
FlowDriver.ResetConvergence()
199+
FlowDriver.Preprocess(0)
200+
FlowDriver.Run()
201+
FlowDriver.Postprocess()
202+
FlowDriver.Update()
203+
stopCalc = FlowDriver.Monitor(0)
204+
```
205+
206+
We need to recover the flow loads and apply them to the structural simulation in order to run the primal iteration for the recording,
207+
208+
```
209+
flow_loads=[]
210+
for iVertex in range(nVertex_Marker_Flow):
211+
vertexLoad = FlowDriver.GetFlowLoad(FlowMarkerID, iVertex)
212+
flow_loads.append(vertexLoad)
213+
214+
FEADriver.SetFEA_Loads(FEAMarkerID, 0, flow_loads[1][0], flow_loads[1][1], flow_loads[1][2])
215+
FEADriver.SetFEA_Loads(FEAMarkerID, 1, flow_loads[0][0], flow_loads[0][1], flow_loads[0][2])
216+
for iVertex in range(2, nVertex_Marker_FEA):
217+
FEADriver.SetFEA_Loads(FEAMarkerID, iVertex, flow_loads[iVertex][0], flow_loads[iVertex][1], flow_loads[iVertex][2])
218+
```
219+
220+
and also, we need to extract the cross dependency on the mesh displacements
221+
222+
```
223+
flow_sensitivities=[]
224+
for iVertex in range(nVertex_Marker_Flow):
225+
sensX, sensY, sensZ = FlowDriver.GetMeshDisp_Sensitivity(FlowMarkerID, iVertex)
226+
flow_sensitivities.append([sensX, sensY, sensZ])
227+
```
228+
229+
that will be used as a source term into the structural domain
230+
231+
```
232+
FEADriver.SetSourceTerm_DispAdjoint(FEAMarkerID,0,flow_sensitivities[1][0],flow_sensitivities[1][1],flow_sensitivities[1][2])
233+
FEADriver.SetSourceTerm_DispAdjoint(FEAMarkerID,1,flow_sensitivities[0][0],flow_sensitivities[0][1],flow_sensitivities[0][2])
234+
for iVertex in range(nVertex_Marker_FEA):
235+
FEADriver.SetSourceTerm_DispAdjoint(FEAMarkerID,iVertex,flow_sensitivities[iVertex][0],flow_sensitivities[iVertex][1],flow_sensitivities[iVertex][2])
236+
```
237+
238+
Next, the structural adjoint simulation is run with
239+
240+
```
241+
FEADriver.ResetConvergence()
242+
FEADriver.Preprocess(0)
243+
FEADriver.Run()
244+
FEADriver.Postprocess()
245+
FEADriver.Update()
246+
stopCalc = FEADriver.Monitor(0)
247+
```
248+
249+
and the crossed sensitivities with respect to the flow load are retrieved using
250+
251+
```
252+
fea_sensitivities=[]
253+
for iVertex in range(nVertex_Marker_FEA):
254+
sensX, sensY, sensZ = FEADriver.GetFlowLoad_Sensitivity(FEAMarkerID, iVertex)
255+
fea_sensitivities.append([sensX, sensY, sensZ])
256+
```
257+
258+
Finally, these boundary displacements are imposed to the flow domain in the next iteration. Once the loop is completed, it only remains to write the solution of each domain to file using
259+
260+
```
261+
FlowDriver.Output(0)
262+
FEADriver.Output(0)
263+
```
264+
265+
### Running SU2
266+
267+
Follow the links provided to download the [python script](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/run_fsi_adjoint.py) and the two sub-config files for the [flow](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/config_channel_ad.cfg) and [structural](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/config_cantilever_ad.cfg) subproblems.
268+
269+
Also, you will need the two mesh files for the [flow domain](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/mesh_channel.su2) and the [cantilever](https://github.com/su2code/Tutorials/tree/feature_python_fsi/multiphysics/python_fsi/mesh_cantilever.su2).
270+
271+
Execute the code using Python
272+
273+
```
274+
$ python run_fsi_adjoint.py
275+
```
276+
277+
The convergence history of each individual domain will be printed to screen, starting with the flow adjoint simulation
278+
279+
```
280+
-------------------------------------------------------------------------
281+
Direct iteration to store the primal computational graph.
282+
Compute residuals to check the convergence of the direct problem.
283+
log10[U(0)]: -16.7835, log10[U(1)]: -16.1388, log10[U(2)]: -17.2225.
284+
log10[U(3)]: -32.
285+
-------------------------------------------------------------------------
286+
287+
+---------------------------------------------------+
288+
| Inner_Iter| rms[A_P]| rms[A_U]| Sens_Geo|
289+
+---------------------------------------------------+
290+
| 0| -2.613161| -1.287093| 0.0000e+00|
291+
| 10| -3.599787| -2.867742| 0.0000e+00|
292+
| 20| -4.291071| -3.771486| 0.0000e+00|
293+
| 30| -5.000513| -4.533287| 0.0000e+00|
294+
| 40| -5.721670| -5.350516| 0.0000e+00|
295+
| 50| -6.457920| -5.924181| 0.0000e+00|
296+
| 60| -7.217678| -6.566948| 0.0000e+00|
297+
| 70| -7.978391| -7.212931| 0.0000e+00|
298+
| 80| -8.738623| -7.908635| 0.0000e+00|
299+
| 90| -9.463444| -8.636672| 0.0000e+00|
300+
| 100| -10.163400| -9.400965| 0.0000e+00|
301+
| 110| -10.847301| -10.206336| 0.0000e+00|
302+
| 120| -11.539210| -11.037515| 0.0000e+00|
303+
| 130| -12.246150| -11.854391| 0.0000e+00|
304+
| 137| -12.812627| -12.143415| 0.0000e+00|
305+
306+
Recording the computational graph with respect to the secondary variables.
307+
```
308+
309+
followed by the structural domain
310+
311+
```
312+
-------------------------------------------------------------------------
313+
Direct iteration to store the primal computational graph.
314+
Compute residuals to check the convergence of the direct problem.
315+
UTOL-A: -14.5125, RTOL-A: -11.6402, ETOL-A: -28.601.
316+
-------------------------------------------------------------------------
317+
318+
+----------------------------------------------------------------+
319+
| Inner_Iter| rms[Ux_adj]| rms[Uy_adj]| Sens[E]| Sens[Nu]|
320+
+----------------------------------------------------------------+
321+
| 0| 1.608195| 2.017211| 1.4163e-05| 6.1467e-01|
322+
| 1| -7.838855| -7.470968| 1.4163e-05| 6.1467e-01|
323+
| 2| -7.472256| -7.061584| 1.4163e-05| 6.1467e-01|
324+
```
325+
326+
We observe that, in both cases, the direct iteration recovers converged residuals for the primal problem.
327+
328+
After 15 iterations, both the flow and structural adjoint fields are successfully converged,
329+
330+
```
331+
+---------------------------------------------------+
332+
| Inner_Iter| rms[A_P]| rms[A_U]| Sens_Geo|
333+
+---------------------------------------------------+
334+
| 0| -11.703018| -10.364366| 5.3425e+02|
335+
| 10| -12.847050| -12.112670| 5.3425e+02|
336+
337+
+----------------------------------------------------------------+
338+
| Inner_Iter| rms[Ux_adj]| rms[Uy_adj]| Sens[E]| Sens[Nu]|
339+
+----------------------------------------------------------------+
340+
| 0| -8.235703| -8.050333| 1.1606e-05| 5.0401e-01|
341+
| 1| -7.792824| -7.398341| 1.1606e-05| 5.0401e-01|
342+
| 2| -7.994688| -7.626300| 1.1606e-05| 5.0401e-01|
343+
```
344+
345+
with a gradient of the drag with respect to the Young's modulus, E, of 1.1606E-05. The convergence of the Young's modulus sensitivity is as follows:
346+
347+
![FSI Results2](../multiphysics/images/fsipyadj1.png)
348+
349+
#### Sensitivity verification
350+
351+
In order to verify the Young's modulus sensitivity, we access the file ```history_ad_1.csv```, which stores the value with 10 significant figures. We obtain
352+
353+
```
354+
Sens[E] = 1.160565822e-05
355+
```
356+
357+
Now, we run central differences to test the accuracy of this value. We compute the drag coefficient for the incremented and decremented value of the Young's modulus
358+
359+
| Young's Modulus $$E$$ | Drag coefficient $$C_D$$ |
360+
|---|---|
361+
| 4.995E+04 | 3.122565491 |
362+
| 5.000E+04 | 3.123146337 |
363+
| 5.005E+04 | 3.123726058 |
364+
365+
which yields a $$\mathrm{d}C_D / \mathrm{d} E$$ computed with central differences of 1.160567E-05, which has an excellent agreement with the value computed via the adjoint.
366+
367+
### Attribution
368+
369+
If you are using this content for your research, please kindly cite the following reference in your derived works:
370+
371+
Sanchez, R. _et al._ (2018), [Coupled Adjoint-Based Sensitivities in Large-Displacement Fluid-Structure Interaction using Algorithmic Differentiation](https://spiral.imperial.ac.uk/handle/10044/1/51023), _Int J Numer Meth Engng, Vol 111, Issue 7, pp 1081-1107_. DOI: [10.1002/nme.5700](https://doi.org/10.1002/nme.5700)
372+
373+
<dl>
374+
This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>
375+
<br />
376+
<a rel="license" href="http://creativecommons.org/licenses/by/4.0/"><img alt="Creative Commons Licence" style="border-width:0" src="https://i.creativecommons.org/l/by/4.0/88x31.png" /></a>
377+
</dl>
378+

0 commit comments

Comments
 (0)