Unsteady aerodynamic optimization problem

[wndzjl]

I'm working on the NACA0012 airfoil with a manageable drag at 15 degrees of attack. Optimization for increased lift. But the optimization process kept reporting errors. Including but not limited to:
① Unsolvable -1 problem,for example:
---- Solving unsteady adjoint for CMZ. t = 3.000000 ----
Calculating dFdW using reverse-mode AD
Using pre-computed KSP PC mat for 3.000000
Solving Linear Equation... 18343 s
Main iteration 0 KSP Residual norm 3.946590696666e+03 18343 s
Main iteration 0 KSP Residual norm 3.946590696666e+03 18343 s
Total iterations 0
Solving Linear Equation... Completed! 18343 s
Residual tolerance not satisfied, solution failed!
Computing total derivatives of CMZ wrt shape
Calculating dFdXv using reverse-mode AD
Calculating [dRdXv]^T * Psi using reverse-mode AD
Computing total derivatives of CMZ wrt aoa
Calculating [dRdAOA]^T * Psi using reverse-mode AD
Calculating [dRdWOld]^T * Psi using reverse-mode AD with level 1
Calculating [dRdWOld]^T * Psi using reverse-mode AD with level 2
Objective Function Sensitivity:
{'DVCon1_volume_constraint_0': {'shape': array([[-0.3188065019420805, -0.3188065019420807, 0.3188065019420806,
..., -0.1238864928490999, 0.1238864928490997,
0.1238864928490994]])}, 'DVCon1_thickness_constraints_0': {'shape': array([[-6.785371384545954 , -0.006792163548094 , 6.785371384546011 ,
..., 0. , 0. ,
0. ],
[-3.3800998287292168, -0.0033834833120413, 3.3800998287292168,
..., 0. , 0. ,
0. ],
[-1.0160393052433474, -0.001017056361605 , 1.0160393052433478,
..., 0. , 0. ,
0. ],
...,
[ 0. , 0. , 0. ,
..., -0.9853716075499095, 0.0009863579655154,
0.9853716075499095],
[ 0. , 0. , 0. ,
..., -3.327114922561015 , 0.0033304453679276,
3.3271149225610177],
[ 0. , 0. , 0. ,
..., -6.716703731000905 , 0.006723427158193 ,
6.71670373100045 ]])}, 'CD': {'shape': array([59.09265063825851 , 61.19576386963102 , 15.698543182754833, ...,
7.379131279394333, 34.52864586018066 , 14.96371529016889 ]), 'aoa': array([-1.4486275768011652])}, 'CL': {'shape': array([2779.457280527302 , 2878.4081822285307, 738.472146470823 , ...,
347.0665353279569, 1624.2177678366884, 703.8443663053816]), 'aoa': array([-68.15762891618162])}, 'CMZ': {'shape': array([-27.497055380702907 , -28.474922952335653 , -7.303340541359451 ,
..., -3.4346530933256485, -16.065235258702742 ,
-6.963451363528834 ]), 'aoa': array([0.6737994965823948])}, 'fail': True}
Adjoint Runtime: 12.0741 s
/home/dafoamuser/dafoam/packages/miniconda3/lib/python3.8/site-packages/pyoptsparse/pySNOPT/pySNOPT.py:495: ComplexWarning: Casting complex values to real discards the imaginary part
hs, xs, pi, rc, inform, mincw, miniw, minrw, nS, ninf, sinf, obj = snopt.snkerc(
Traceback (most recent call last):
File "runScript_v2.py", line 260, in
sol = opt(optProb, sens=optFuncs.calcObjFuncSens, storeHistory=histFile)
File "/home/dafoamuser/dafoam/packages/miniconda3/lib/python3.8/site-packages/pyoptsparse/pySNOPT/pySNOPT.py", line 524, in call
sol_inform["text"] = self.informs[inform]
KeyError: -1
Finalising parallel run
Finalising parallel run

② Gradient divergence problem: The same running code, add an FFD point in the x direction, or encrypt the grid a little, there is a gradient calculation problem. I tried to change "pcFillLevel": 1 to 2 or even 3, but of course I don't know what the upper limit is. And the time step is encrypted. It's all not quite right.
Here is the gradient before the change:
"Optimization_Iteration_001":
{
"CD":
{
"aoa": [ 3.631911043273110e-02 ],
"shape": [ 1.146847873014347e+00, 1.147068502241482e+00, 3.900431154067007e-01, 3.905247864378928e-01, -2.810339625258564e-02, -2.805524484318332e-02, -6.708802728112315e-01, -6.703186697641548e-01, 1.822455757345192e-02, 1.794116130871638e-02, -3.925997648753567e-01, -3.938004630533021e-01, 2.354189400231279e-01, 2.346852036960353e-01, 1.956921489772387e-02, 1.748167234558832e-02, -2.352356102830359e-01, -2.357947779020004e-01, -4.809955587673668e-01, -4.820212683928558e-01 ]
},
"CL":
{
"aoa": [ -1.368855429717181e-01 ],
"shape": [ -1.041966873608299e+01, -1.042859475811789e+01, -2.678775406386524e+00, -2.695648919194490e+00, 8.097231047151267e-01, 8.046943777414435e-01, 6.244354965956288e+00, 6.217138681941133e+00, 7.001332643883140e-01, 6.966902017413777e-01, 3.762593056302457e+00, 3.748356755425958e+00, 7.588840263699398e-01, 7.537621263061530e-01, 9.599373588052813e-01, 9.458432414730797e-01, 1.899495234664280e-01, 1.828042026119045e-01, -2.709149390797907e-01, -2.812621283832193e-01 ]
},
"CMZ":
{
"aoa": [ -1.885842074905516e-02 ],
"shape": [ -3.674070510526958e-01, -3.674873993680131e-01, -8.808975246393418e-02, -8.778858585820232e-02, 4.796934097658802e-02, 4.812967041650476e-02, 3.018032861861055e-01, 3.024528008755866e-01, -1.209706399169061e-01, -1.204810946392721e-01, 1.894166312662931e-01, 1.913689224328309e-01, -6.468927526942483e-01, -6.459116908419558e-01, -1.733880742994205e-01, -1.704713731691231e-01, 8.928265726047417e-02, 8.973439880976847e-02, 7.056085770910947e-01, 7.067027136236040e-01 ]
}
},

The following is the gradient of the same iteration step after modification:
"Optimization_Iteration_001":
{
"CD":
{
"aoa": [ -2.641444438970758e+02 ],
"shape": [ -5.954049781959720e+03, -5.955545067645521e+03, -3.704027927802165e+03, -3.711359363861115e+03, -3.757427539027314e+02, -3.767795152323452e+02, 1.396735981333199e+03, 1.384301229572612e+03, 5.492533218738623e+02, 5.509717002039683e+02, 5.774524331711179e+03, 5.783948076917006e+03, 1.726971576369909e+02, 1.710689184339229e+02, 1.490943258412250e+03, 1.490837920382954e+03, -1.744425641379867e+02, -1.809253521041125e+02, -2.214876079005413e+01, -3.230845164116998e+01, 3.934509259939206e+02, 3.872099631777002e+02, 4.742213440211185e+02, 4.671654094067212e+02 ]
},
"CL":
{
"aoa": [ -1.242929701025411e+04 ],
"shape": [ -2.801987326371587e+05, -2.802691087996546e+05, -1.743004943361459e+05, -1.746454912740329e+05, -1.767543567579631e+04, -1.772424801116990e+04, 6.574168028187593e+04, 6.515642860084804e+04, 2.584514335225823e+04, 2.592596289363368e+04, 2.717204604011194e+05, 2.721636301079290e+05, 8.123901637884427e+03, 8.047317069910652e+03, 7.015360989561443e+04, 7.014863278064877e+04, -8.216678735332131e+03, -8.521619835263229e+03, -1.047801911893095e+03, -1.525719298900452e+03, 1.852899161725804e+04, 1.823546852084492e+04, 2.233298241156832e+04, 2.200112094403373e+04 ]
},
"CMZ":
{
"aoa": [ 1.228521291945287e+02 ],
"shape": [ 2.769152664797490e+03, 2.769847949911460e+03, 1.722753566080214e+03, 1.726163177695115e+03, 1.748311747742761e+02, 1.753130063958440e+02, -6.495356706888278e+02, -6.437545277175977e+02, -2.554640585614549e+02, -2.562636218754770e+02, -2.685622418373853e+03, -2.690008340287385e+03, -8.047126910721288e+01, -7.971339531069124e+01, -6.934662317065855e+02, -6.934170141910200e+02, 8.065716506476200e+01, 8.367365442116946e+01, 1.000900788765625e+01, 1.473612505523211e+01, -1.827353562311075e+02, -1.798307267632012e+02, -2.201319675936968e+02, -2.168481855422431e+02 ]
}

[friedenhe]

Can you add "useMGSO": True, to adjEqnOption in runScript_v2.py to see if this improves the adjoint convergence?

[wndzjl]

Yeah, I tried that, too. But the results are still not good. Previous simulation I set "useNonZeroInitGuess": True, but as you can see below, the "Main iteration 0 KSP Residual norm"Getting bigger and bigger. In theory, the non-zero assumption is used, and the exact solution of the previous step is the initial solution of the next step. The difference between one step and the other can't be that far. It's basically an order of magnitude growth rate. After it is controlled by relative convergence, ksp is still large after iteration, and the gradient diverges.
---- Solving unsteady adjoint for CL. t = 2.892000 ----
Calculating dFdW using reverse-mode AD
Updating dRdWTPC mat value using OF fvMatrix
Solving Linear Equation... 37700 s
Main iteration 0 KSP Residual norm 2.553806483028e+99 37701 s
Main iteration 100 KSP Residual norm 4.314409305259e+86 37705 s
Main iteration 200 KSP Residual norm 1.904602111316e+55 37710 s
Main iteration 300 KSP Residual norm 1.247276594790e+17 37716 s
Main iteration 343 KSP Residual norm 2.057444609526e-01 37718 s
Total iterations 343
Solving Linear Equation... Completed! 37718 s
Residual tolerance satisfied, solution finished!
Computing total derivatives of CL wrt shape
Calculating dFdXv using reverse-mode AD
Calculating [dRdXv]^T * Psi using reverse-mode AD
Computing total derivatives of CL wrt aoa
Calculating [dRdAOA]^T * Psi using reverse-mode AD
Calculating [dRdWOld]^T * Psi using reverse-mode AD with level 1
Calculating [dRdWOld]^T * Psi using reverse-mode AD with level 2
---- Solving unsteady adjoint for CL. t = 2.891000 ----
Calculating dFdW using reverse-mode AD
Updating dRdWTPC mat value using OF fvMatrix
Solving Linear Equation... 37721 s
Main iteration 0 KSP Residual norm 2.542383510874e+100 37721 s
Main iteration 100 KSP Residual norm 8.736555995282e+86 37726 s
Main iteration 200 KSP Residual norm 2.019323637527e+55 37731 s
Main iteration 300 KSP Residual norm 7.464065924228e+16 37737 s
Main iteration 340 KSP Residual norm 1.606598197069e+00 37739 s
Total iterations 340
Solving Linear Equation... Completed! 37739 s
Residual tolerance satisfied, solution finished!
Computing total derivatives of CL wrt shape
Calculating dFdXv using reverse-mode AD
Calculating [dRdXv]^T * Psi using reverse-mode AD
Computing total derivatives of CL wrt aoa
Calculating [dRdAOA]^T * Psi using reverse-mode AD
Calculating [dRdWOld]^T * Psi using reverse-mode AD with level 1
Calculating [dRdWOld]^T * Psi using reverse-mode AD with level 2
---- Solving unsteady adjoint for CL. t = 2.890000 ----
Calculating dFdW using reverse-mode AD
Updating dRdWTPC mat value using OF fvMatrix
Solving Linear Equation... 37741 s
Main iteration 0 KSP Residual norm 2.521157307970e+101 37741 s
Main iteration 100 KSP Residual norm 2.702167467956e+88 37746 s
Main iteration 200 KSP Residual norm 6.902745573252e+56 37751 s
Main iteration 300 KSP Residual norm 2.456498008163e+18 37757 s
Main iteration 341 KSP Residual norm 1.655384038143e+01 37759 s
Total iterations 341
Solving Linear Equation... Completed! 37759 s
Residual tolerance satisfied, solution finished!
Computing total derivatives of CL wrt shape
Calculating dFdXv using reverse-mode AD
Calculating [dRdXv]^T * Psi using reverse-mode AD
Computing total derivatives of CL wrt aoa
Calculating [dRdAOA]^T * Psi using reverse-mode AD
Calculating [dRdWOld]^T * Psi using reverse-mode AD with level 1
Calculating [dRdWOld]^T * Psi using reverse-mode AD with level 2
logOpt.txt

[wndzjl]

In the same case, I changed the time step from 0.001 to 0.002, at least in the same number of reverse derivations, without gradient divergence or 0 ksp being too large. This is a little incomprehensible, why the time step is larger, but the effect is better. Thank you！

[friedenhe]

Yes, you may need to play with the time step. Usually, reducing the time step should make the case work better. I am surprised to see that increasing the time step works better.....

[wndzjl]

Thank you for your reply! What is the relationship between the initial residual "0 KSP residual" and the total physical time for the same situation and time step? For example, when dt=0.001, is the convergence effect of simulating physics time 0-2 seconds better than 0-1 seconds, even though the former takes more time? I think I have tried such a case before and would like to see if you have any experience in this area?

[friedenhe]

initial residual "0 KSP residual" is for adjoint computation and has nothing to do with the physical time in flow simulations. Not sure what you meant in your second question...

[wndzjl]

I'm sorry I didn't make myself clear. Does the length of simulation time affect convergence in the same case?

[friedenhe]

No, the length of simulation time does NOT impact the convergence.

[wndzjl]

Dear Professor He, I am sorry to ask you again about the initial residual 0 KSP residual. All I've seen so far is that the adjoint equations converge, but the calculated gradients are unusually large. The direct reason is that in the process of reverse derivation, at a certain time, the initial residual 0 KSP residual increases by orders of magnitude in the process of downward calculation. This usually does not occur in the first gradient calculation, that is, in the process of solving the gradient after the deformation. I have tried to reduce the time step or modify the tolerance in the fvSolution, and I was surprised to find that reducing the tolerance, for example, from the original 1e-8 to 1e-10, would accelerate the divergence process. I would like to ask you again what other parameters this initial residual 0 KSP residual is related to and how I can improve it. I look forward to your reply and wish you a happy New Year！

[friedenhe]

We are 100% sure why this happened to you. Have you tried other cases to see if this happens for all cases? Also, can you try to use a loose tolerance, like 1e-6, helps with the convergence?

[wndzjl]

No,I've already tried the 1e-6.This is not a very good case, but there are successful cases. I just don't know why this case doesn't work