Extracting Residual Jacobians to Conduct a Custom Adjoint Optimization

[Bot-Enigma-0]

Hi,

I am working on an adjoint optimisation project using SU2 as the simulation tool and JAX for automatic differentiation. My goal is to optimise a custom objective function related to the boundary layer, which requires computations that (to my knowledge) cannot be done in SU2. I have successfully used JAX to compute the derivative of the objective function with respect to the flow field (∂J/∂U). To form the adjoint problem, I need the derivatives of the residual with respect to both the flow field and the design variables (∂R/∂U and ∂R/∂x).

From my understanding, the derivatives, particularly ∂R/∂U, seem to be stored in the source files CSysSolve_b.cpp and CSysMatrix.cpp. I am unsure about the location of ∂R/∂x. I want to know how to extract these Jacobians? Specifically:

1. Is there an existing method to access or extract the Jacobians from these files, or would it require modifications to the source code?
2. If modifications are needed, should they be made within the CSysMatrix.cpp file, or are there additional changes required in other files, such as the solver-specific ones (e.g., CAdjNSSolver.cpp)?
3. My knowledge in C++ is limited, so any suggestions are welcome!

Thanks
RK

[pcarruscag]

The Jacobian in CSysMatrix is an approximation, it won't give you correct derivatives.
Products with the true Jacobian are computed with AD, we do not store the true Jacobian anywhere.

Why don't you modify SU2 to compute the objective function you want?
What is the objective function btw?

[Bot-Enigma-0]

Hi,

Sorry for the late reply. I editted the code a bit, I think this should work (https://github.com/Bot-Enigma-0/SU2-Inject-Custom-Grad). Running this code gives finite gradient values for the design variables, however they are quite large. I'm not sure if this is because the objective function is relating to the boundary layer, where the gradients are expected to be large. Deforming the grid with these results in the entire domain to be affected, the shape is quite a mess; the gradient norm in this case is 1.930993856 (left). Damping these gradients with resepct to the boundary layer thickness gives something more reasonable (and expected). Here the gradient norm is 0.000420048 (right); I have zoomed into one of the deformations:

However a few iterations of doing this causes the solution to diverge; the deformations get bigger, even when damped. I think this is because the line search between the 'major' design iterations is not being done. Is this meant to be done outside SU2? I understand that SU2 uses scipy for the optimisation process when shape_optimisation.py is called. Will I need to write my own optimisation module for my case? If so will writing a python wrapper help? There were some points in #923 to use python wrappers for the optimisation process, not sure if this has been worked on since then. I think this would also help with automating the process. Currently I am manually executing each script individually.

To validate the changes to the code, I have written a test case for an inviscid wedge to optimise for lift and drag (https://github.com/Bot-Enigma-0/Testcase_for_custom_gradients). The results show good agreement with the deformation of the surface (especially for the drag). The adjoint fields are completely wrong, even though residuals for both dropped to -14, -16, -17 and -19 orders of magnitude for rho, rhoU, rhoV and rho*E. The sensitivity fields look similar but closer inspection of each sensitivity shows is a big difference at some places. If the value range for the 'custom' method is set to the range of the SU2 given adjoint fields, they appear to share similar contours but the custom method in general has greater values. I'm not sure if there is something I am missing in my formulation of the gradients, which are causing such large differences in the adjoint fields. For all of the objective functions I have written, there is a large value at the inlets and outlets in the adjoint fields.

Here are some brief results for the validation:
Flow Variables:

	SU2	custom	% diff
Lift	-1.568438	-1.543804	1.5706
Drag	0.387356	0.382944	1.1390

Deformation Gradient norm:

	SU2	custom	% diff
Lift	0.05166527	0.04010482	-22.37566
Drag	0.01262895	0.01149017	-9.017218

The grid deformation for lift and drag:

I have many more results comparing adjoint fields, sensitivities and FFD gradients here (https://github.com/Bot-Enigma-0/Testcase_for_custom_gradients/tree/main/compare), posting it here will make the post quite long.

Thanks,
RK

[pcarruscag]

You can use https://github.com/su2code/FADO, see the included examples with SU2.

[Bot-Enigma-0]

Hi,

I've been playing around with FADO, writing scripts for some tutorials. I kind of get the gist of it, its pretty much exactly what I had in mind. Example 4 is for a CFD optimisation problem. In this all of the optimisation parameters (objecive function and constraints) are manually called from this history output. If the objective function is 'manually' set inside the FADO script, I guess I can write my objective function in this script, or have it called.

For the lift and drag testcase I mentioned in my last post, the adjoint fields for the custom/coupled objective functions seem quite 'noisy', but have some regions of resemblance. The gradient norm, sensitivities and deformation are reasonably close. Are there any additional terms I should inject into the source code? Currently I am finding and injecting ∂J/∂U and ∂J/∂X in the source code. The adjoint converge for the both is quite good - 14, -16, -17 and -19 orders of magnitude for rho, rhoU, rhoV and rho*E; the primal convergence was along these lines as well. Here is a comparison of the adjoint fields for drag (the first is set to the default range, and the second scaled to the range of the SU2 given field):

More of the results are here: https://github.com/Bot-Enigma-0/Testcase_for_custom_gradients/tree/main/compare

Thanks,
RK

[pcarruscag]

The definite proof is verifying the gradients you obtain from SU2, since you have very good convergence of the primal you should be able to use finite differences.

[Bot-Enigma-0]

Hi, I've done the validation of 20 design variables for step sizes ranging from e-1 to e-10 for both the lift and drag cases. there seems to be an acceptable error/difference between the coupled adjoint gradients and the su2 given gradients. Below are the plots comparing the errors between the coupled (finite difference) gradients with the SU2 AD adjoint gradients for drag:

coupled FD against coupled AD:

and SU2 FD gradients against SU2 AD gradients:

The question is, if the FD gradients agree reasonably, the surface sensitivities and sensitivity fields, the AD (dJ/d) gradients, and the deformation given by the coupled method agree quite well (see https://github.com/Bot-Enigma-0/Testcase_for_custom_gradients/tree/main/compare), why would the adjoint field be so noisy and different in the coupled method compared to the SU2 given fields?

Thanks,
RK

[pcarruscag]

I have no idea what you mean by "coupled" or what the variables in those plots are supposed to show.
But if the derivatives agree then it's fine. The adjoint variables are not expected to be smooth with the adjoint equations we solve.

[Bot-Enigma-0]

By coupled I mean the gradients dJ/dU and dJ/d which are found externally and injected into the source code. The plots are showing the errors of each of the 20 design variables wrt the step size.

But if the derivatives agree then it's fine. The adjoint variables are not expected to be smooth with the adjoint equations we solve.

I see, thanks.