[kmandar]

Dear SU2 developers,

I am performing Continuum Sensitivity Analysis with SU2 to output shape sensitivities. This gives me direct (forward) results. For example, with the NACA0012 quick start problem, I can get derivatives of pressure at all ponits on the airfoil with respect to change in shape of the NACA0012 airfoil (with Hicks-Henne design variable).

Is there a way to get pressure derivatives on the airfoil so that I can verify the direct results?

Can Algorithmic Differentiation (AD) implemented in version 4.2 be of any use to me?

Using the SU2 adjoint, I can compare the derivative of Lift of Drag (i.e. integral of pressure over the airfoil), but not the derivative of pressure on the airfoil itself. Do you have any comments on this?

Thanks for your help,
Mandar

[talbring]

Hi kmandar,

what do you mean with pressure derivatives and direct results, did you implement CSA in SU2 ?

To get the gradient of the pressure at each point on the surface with respect to a design variable you have to use the direct differentiation mode and add some output in the code itself.

First of all see here for Direct Differentiation mode (Forward mode of AD) compilation:
http://www.cfd-online.com/Forums/su2...ifference.html

When running with direct diff. mode enabled (DIRECT_DIFF=DESIGN_VARIABLES in the config file and using SU2_CFD_DIRECTDIFF as the executable) you can get the gradient of any variable inside the code with respect to a design variable with SU2_TYPE::GetDerivative(..). The particular design variable is specified with the usual DV_KIND, DV_MARKER, DV_PARAM options and DV_VALUE must be set to a non-zero value.

[kmandar]

Hi Tim,

I am using SU2 as a black-box to get the sensitivities using CSA. For this however, I do need the Jacobian matrix that is computed in SU2 during Euler implicit time marching. So, I have made minor changes in SU2 to output this Jacobian matrix. The process to get sensitivities is: (a) I run SU2_CFD to get the flow analysis results and also output the Jacobian matrix, (b) I solve a linear system of equations external to SU2 to get the direct sensitivities.

I appreciate your reply regarding the direct differentiation mode. I am trying that out.

Why is it needed for the design variable to have a non-zero value (non-zero DV_VALUE)? I would assume that to get a derivative, only the design variable definition should be enough (eg. Hicks-Henne bump at mid-chord on the upper surface).

Thanks,
Mandar

[talbring]

Originally Posted by kmandar

Why is it needed for the design variable to have a non-zero value (non-zero DV_VALUE)? I would assume that to get a derivative, only the design variable definition should be enough (eg. Hicks-Henne bump at mid-chord on the upper surface).

It is implemented like this because of the python script. The script essentially loops through all specified design variables and sets DV_VALUE for one design variable to 1 and for all others to 0. If you have more than one DV_VALUE set to one, you get the sum of the gradients.

[kmandar]

Hi talbring,

That answers my question.

Thanks,
Mandar

[kmandar]

Hi talbring,

I have two separate questions, but both related to the forward differentiation mode.

(1) Use of SU2_CFD_DIRECTDIFF executable

Is there an test case for the use of SU2_CFD_DIRECTDIFF?

After I ran SU2_CFD_DIRECTDIFF, I get the derivatives of all the force coefficients in the history file. However, I am more interested in getting the derivatives of the conservative and primitive flow variables using the direct differentiation mode.

Originally Posted by talbring

When running with direct diff. mode enabled (DIRECT_DIFF=DESIGN_VARIABLES in the config file and using SU2_CFD_DIRECTDIFF as the executable) you can get the gradient of any variable inside the code with respect to a design variable with SU2_TYPE::GetDerivative(..). The particular design variable is specified with the usual DV_KIND, DV_MARKER, DV_PARAM options and DV_VALUE must be set to a non-zero value.

Can you please give some details about how I can get the gradient of any variable (for example pressure or density or Conservative_2) with respect to a design variable? I could not find the function GetDerivative or the class SU2_TYPE on the SU2 doxygen page. Where should I find this?

Also, since these derivatives would be at each point in the domain, I believe these cannot be output in the history file, but instead in a file similar to flow.dat.

(2) Use of the python script

I also tried to run the direct_differentiation.py python scrip as you mentioned in this thread. I run it as:

direct_differentiation.py -n 16 -f inv_NACA0012_DIRECTDIFF_h1.cfg > DIRECT_DIFF_PyScript_log.txt

but I get the following error:

KeyError: 'Config parameter not found: DEFINITION_DV'

Am I missing something here? Also is it correct to run this with multiple processors using

-n 16

as I have done above?

Thanks a lot,
Mandar

[kmandar]

I got the answer to my second question -- I had not specified my Hicks-Henne design varible in the CFG file as
DEFINITION_DV= ( 1, 1.0 | airfoil | 1, 0.50 )
Once this is done, the python script executes without any issues.

However, I do not understand why we need to define the design variable at two places in the CFG file:

- once in the section "DESIGN VARIABLE PARAMETERS" with
DV_KIND= HICKS_HENNE
DV_MARKER= ( airfoil )
DV_PARAM= ( 1, 0.5 )
DV_VALUE= 0.01

and again

- in the section "OPTIMAL SHAPE DESIGN DEFINITION" with
DEFINITION_DV= ( 1, 1.0 | airfoil | 1, 0.50 )

Since both these are same, won't specifying the design variable at either of the places be enough?

Thanks,
Mandar

[talbring]

Originally Posted by kmandar

Hi talbring,

I have two separate questions, but both related to the forward differentiation mode.

(1) Use of SU2_CFD_DIRECTDIFF executable

Is there an test case for the use of SU2_CFD_DIRECTDIFF?

After I ran SU2_CFD_DIRECTDIFF, I get the derivatives of all the force coefficients in the history file. However, I am more interested in getting the derivatives of the conservative and primitive flow variables using the direct differentiation mode.



Can you please give some details about how I can get the gradient of any variable (for example pressure or density or Conservative_2) with respect to a design variable? I could not find the function GetDerivative or the class SU2_TYPE on the SU2 doxygen page. Where should I find this?

Also, since these derivatives would be at each point in the domain, I believe these cannot be output in the history file, but instead in a file similar to flow.dat.

Yes you are right, they would need to be in flow.dat.
This is a feature I definitely want to implement in the future. But I haven't had the time to do this.

Also note that the doxygen page is based on a quite old version of SU2. It is better to look at the actual source code or to create the doxygen documentation from the Documentation repository. Then you also find the GetDerivative routine that you can use on the primitive or conservative variables.

Originally Posted by kmandar

However, I do not understand why we need to define the design variable at two places in the CFG file:

- once in the section "DESIGN VARIABLE PARAMETERS" with
DV_KIND= HICKS_HENNE
DV_MARKER= ( airfoil )
DV_PARAM= ( 1, 0.5 )
DV_VALUE= 0.01

and again

- in the section "OPTIMAL SHAPE DESIGN DEFINITION" with
DEFINITION_DV= ( 1, 1.0 | airfoil | 1, 0.50 )

Since both these are same, won't specifying the design variable at either of the places be enough?

DV_KIND, DV_MARKER, DV_PARAM and DV_VALUE are the options for the binaries (e.g. for SU2_DEF and SU2_DOT/SU2_DOT_AD) and DEFINITION_DV is used to give the python scripts the information of the design variables. Inside the python framework this information will be unpacked to set the DV_* options for the individual calls of the binaries. So the actual values of the DV_* options are not important if you run the python scripts because they will be overwritten anyway.

[kmandar]

Hi Tim,

Originally Posted by talbring

Yes you are right, they would need to be in flow.dat.
This is a feature I definitely want to implement in the future. But I haven't had the time to do this.

Ok, I will stay tuned for updates about this!!

I found function calls in solver_adjoint_discrete.cpp such as:

Local_Sens_Press = SU2_TYPE::GetDerivative(Pressure)

;

However, from your comment, it seems that this has to do with pressure at the farfield. Is this correct?

And also some others in output_structure.cpp, but I believe those are just to log the output.

Thanks for the explanation about design variable definition.

Mandar

[kmandar]

Dear SU2 developers,

I am using the SU2 Release 4.1.1 "Cardinal" to calculate the derivative of lift coefficient (CL) for the quick start example (NACA0012) using different methods and the results are as shown below.

I did the analysis for smooth flow over the airfoil (no shock):
MACH_NUMBER= 0.5
AoA= 1.25
FREESTREAM_PRESSURE= 101325.0
FREESTREAM_TEMPERATURE= 288.15

Values of CL derivative with respect to Hicks-Henne bump function on upper surface at mid-chord location i.e DEFINITION_DV= ( 1, 1.0 | airfoil | 1, 0.50 ):

- Discrete Adjoint (run script discrete_adjoint.py): 2.83854
- Discrete Direct (run script direct_differentiation.py): 2.9750073914
- Continuous Adjoint (run SU2_CFD with DIRECT, CONTINUOUS_ADJOINT, then run SU2_DOT): 2.86376

(Q1) Why should the Discrete Adjoint and Discrete Direct results differ?
(Q2) Which one of these should I trust the most?


I was assuming that if I refine the mesh enough, all the above results would converge to the same (or at least similar) value. So, I refined the mesh uniformly, i.e each triangle is split into 4 triangles by joining the midpoints of the three sides. The original mesh has ~10k triangles. The finesh mesh that I have (obtained after 4 such refinements) has ~10k*(4^4) = 2600k i.e. ~2.5 million triangles. For this mesh, the three values are:

- Discrete Adjoint (run script discrete_adjoint.py): 2.83487
- Discrete Direct (run script direct_differentiation.py): 2.971274018
- Continuous Adjoint (run SU2_CFD with DIRECT, CONTINUOUS_ADJOINT, then run SU2_DOT): 2.84072

I have attached a plot which shows results at all the mesh levels.

(Q3) Shouldn't all the results (or at least the Discrete Adjoint and Discrete Direct results) converge to a similar (if not the same) point?

Here are some details from the configuration file, which are same for all the mesh levels:

CONV_NUM_METHOD_FLOW= ROE
SPATIAL_ORDER_FLOW= 2ND_ORDER
TIME_DISCRE_FLOW= EULER_IMPLICIT
MGLEVEL= 3
CONV_CRITERIA= RESIDUAL
RESIDUAL_REDUCTION= 10
RESIDUAL_MINVAL= -10
LINEAR_SOLVER= FGMRES
LINEAR_SOLVER_ERROR= 1E-6
LINEAR_SOLVER_ITER= 5
% For Adjoint method
CONV_NUM_METHOD_ADJFLOW= JST
AD_COEFF_ADJFLOW= ( 0.15, 0.5, 0.02 )
SPATIAL_ORDER_ADJFLOW= 2ND_ORDER
TIME_DISCRE_ADJFLOW= EULER_IMPLICIT

(Q4) What is the expected rate of convergence of shape derivatives (such as the one I am calculating) using each of these methods, with SPATIAL_ORDER_FLOW= 2ND_ORDER?


Thanks for your help,
Mandar

[talbring]

Hi Mandar,

Q1:
I just checked it myself, I get more or less the same numbers for all three methods, even on the coarse mesh. (Discrete Adjoint/ DirectDiff/ Fin. Diff.: 2.97, Cont. Adj.: 3.03)
Are you sure you used OBJECTIVE_FUNCTION=LIFT for the adjoint methods ?

Q2:
In my opinion you can order all four methods (discrete adj. cont. adj, direct diff, finite diff.) like this (of course this is also case dependent):

Accuracy: Direct Diff > Disc. Adj > Cont. Adj > Finite Diff.
Efficiency: Cont. Adj. > Disc. Adj. > Finite Diff. > Direct Diff.

Q3:
Yes, in theory (although there is no proof of that yet) they should converge to the same values for an infinitely fine mesh. However, this is only true if everything is solved exactly, which is never the case (even if the residual is down to machine precision).

Q4:
What convergence do you mean ? The convergence of the derivatives with respect to the accuracy of the flow solution ? Or the convergence of the adjoint solvers itself with respect to the number of iterations?

Tim

[kmandar]

Hi Tim,

I really appreciate your prompt reply.

Q1:
I just checked it myself, I get more or less the same numbers for all three methods, even on the coarse mesh. (Discrete Adjoint/ DirectDiff/ Fin. Diff.: 2.97, Cont. Adj.: 3.03)
Are you sure you used OBJECTIVE_FUNCTION=LIFT for the adjoint methods ?

Yes, I did use OBJECTIVE_FUNCTION=LIFT for the adjoint methods.
Since you are getting the exact same values for Discrete-Adjoint and Discrete-Direct (Direct-diff) methods, I am re-checking my calculations. I understand that it is possible for Continuous adjoint to yield a different result than the discrete methods.

Q2:
In my opinion you can order all four methods (discrete adj. cont. adj, direct diff, finite diff.) like this (of course this is also case dependent):

Accuracy: Direct Diff > Disc. Adj > Cont. Adj > Finite Diff.
Efficiency: Cont. Adj. > Disc. Adj. > Finite Diff. > Direct Diff.

I was expecting Direct Diff = Disc. Adj for accuracy, because they should give the same result but I guess the accuracy may differ due to machine precision issues for a fine mesh.

Q3:
Yes, in theory (although there is no proof of that yet) they should converge to the same values for an infinitely fine mesh. However, this is only true if everything is solved exactly, which is never the case (even if the residual is down to machine precision).

I completely agree with you.

Q4:
What convergence do you mean ? The convergence of the derivatives with respect to the accuracy of the flow solution ? Or the convergence of the adjoint solvers itself with respect to the number of iterations?

I mean the convergence of the derivatives (such as CL derivative) w.r.t. size of the elements (or spatial resolution). Here are some papers: paper-1 (see Figure 4) and paper-2 (see Figures 3 and 4 and Table 4) which show rate (or order) of convergence of shape derivatives.


I have a related question about the Discrete-Adjoint method: Can I output the discrete adjoint vector from the code?

Thanks,
Mandar

[kmandar]

Hi Tim,

I redid the analysis and I have attached my new results. Now the Discrete Adjoint and Direct methods give (almost) same results.

The reason for mismatch was very surprising to me. Consider the following two meshes that differ only in the ordering of nodes on the airfoil. All other things (number of nodes, elements, node locations) are exactly the same.
(a) Node 1 on the airfoil is the first node on the lower surface near the leading edge, and the following nodes go counterclockwise (lower surface, trailing edge, upper surface), ending with Node 200 at the leading edge. This is attached as the mesh file mesh_primary_h1_old.su2.
(b) Node 1 on the airfoil is the first node on the lower surface near the trailing edge, and the following nodes go clockwise (lower surface, leading edge, upper surface), ending with Node 200 at the trailing edge. This is attached as the mesh file mesh_primary_h1.su2. The is same as the ordering for the mesh of the quick-start NACA-0012 example.

I have also attached the input file I used to run these two cases.

With (a), the results for Discrete Adjoint and Direct differ (as given in my previous email), but with (b), they are the same. Do you have an idea of what could be the reason for this?

I zipped the mesh files and CFG file which can be downloaded from this link.

Thanks,
Mandar