Incorrect AP variable function sensitivities for simple 2D shock problem

[garobed1]

Description

I have been using ADFlow to solve a 2D CFD portion of a problem involving a shock impinging on a flat deforming surface, and I am trying to find the sensitivity of the viscous drag on the surface with respect to the shock impingement angle. The shock angle determines the flow properties downstream of the shock (mach, beta, P, T), which I'm trying to get the sensitivies of through the corresponging AeroProblem variables. Beta refers to the angle of attack instead of alpha, since the 2D plane of the mesh is defined on x-z coordinates.

The problem is defined somewhat strangely through the CGNS mesh. The shock originates from the top left corner of the domain. The upstream properties are hard-coded in the CGNS mesh as a supersonic inflow BC on the left surface. The top surface is a far-field BC with properties defined through the AeroProblem, representing the downstream properties of the shock as defined through oblique shock relations. The bottom surface is an inviscid BC up from x= 0 to x=1, and a viscous heat flux wall past that. The right surface is a supersonic outflow BC. As a result, the AeroProblem variables should control only the top surface boundary conditions, allowing for changes in the angle of the shock originating from the top left. The surfaces along the y axis are symmetry planes, and the mesh is 1 cell wide in the y direction.

The derivative of the viscous drag on the bottom surface cdv with respect to mach is correct to 3 places, but the beta, P, and T derivatives are way off. I suspect this is a consequence of the highly irregular mesh definition, but I'm not sure. I would provide the CGNS mesh, but I can't seem to do that in a bug report here.

Steps to reproduce issue

import numpy
from mpi4py import MPI
from baseclasses import *
from adflow import ADFLOW

alpha = 0. #
beta = 7.2833969362749187
mach = 2.6381157549933598
areaRef = 1.0
chordRef = 1.0
T = 254.02071103827234 
P = 4987.6905797938707
probName = 'impinge_mphys'
aeroGridFile = f'./imp_mphys_73_73_25.cgns'

aeroOptions = { #ADflow aero solver options
    # Common Parameters
    'gridFile':aeroGridFile,
    'outputDirectory':'../results/',
    'writeTecplotSurfaceSolution':False,
    'writeSurfaceSolution':True,
    'writeVolumeSolution':True,
    
    # Physics Parameters
    'equationType':'RANS',
    'turbulenceModel':'SA',
    'turbulenceProduction':'vorticity',
    'useft2SA':True,
    'eddyVisInfRatio':3.0,

    # Common Parameters
    "CFL": 1.5,
    "CFLCoarse": 1.25,
    'MGCycle':'sg',
    'nCycles':100000,
    'monitorvariables':["resrho", "resturb"],
    'useNKSolver':True,
    'NKSwitchTol':1e-4,#e-1,
    'NKSubspaceSize':50,
    'NKPCILUFill':3,
    'NKLS':'none',
    'useANKSolver':True,
    'ANKCoupledSwitchTol':1e-3,
    'ANKConstCFLStep':0.4,
    'ANKCFLLimit':1000000000.0,
    "L2Convergence": 1e-12,
    "forcesAsTractions": False,
    
    # Adjoint options
    'adjointL2Convergence': 1e-06,
    # Output
    'volumeVariables':['eddyratio','mach','cp','temp'],
    'surfaceVariables':['yplus','cf','cp','cfx','cfy','cfz'],
    'printIterations':False,
    'printTiming':False,
    'printWarnings':True,
    'setMonitor':False
    }

# Create solver
CFDSolver = ADFLOW(options=aeroOptions)

# Aerodynamic problem description
ap = AeroProblem(
    name=probName,
    mach=mach,
    alpha =alpha,
    beta =beta,
    areaRef = 1.0,
    chordRef = 1.0,
    T = T, 
    P = P, 
    evalFuncs=["cdv"],
)
ap.addDV("mach")    

# Solve and evaluate functions
funcs = {}
funcsens = {}
CFDSolver(ap)
CFDSolver.evalFunctions(ap, funcs)
CFDSolver.evalFunctionsSens(ap, funcsens)


# central difference check, directional derivative
# step size
h = 1e-5

# just do full fd
fd = 0.


funcs2 = {}
ap.mach += h
CFDSolver(ap)
CFDSolver.evalFunctions(ap, funcs2)


funcs3 = {}
ap.mach -= 2*h
CFDSolver(ap)
CFDSolver.evalFunctions(ap, funcs3)

fd = (funcs2['impinge_mphys_cdv'] - funcs3['impinge_mphys_cdv'])/(2*h)

if MPI.COMM_WORLD.rank == 0:
    print(fd)
    print(funcsens['impinge_mphys_cdv'])

Behavior

The only correct gradient that I have found for this configuration is cdv w.r.t. mach. Top number is approximated via central differencing, bottom number is given by evalFunctionsSens

-0.0021169910893520036
OrderedDict([('mach_impinge_mphys', -0.0021172060598884054)])

cdv w.r.t. beta:

4.5211767316088995e-05
OrderedDict([('beta_impinge_mphys', 0.00265107401133413)])

cdv w.r.t. P:

-3.868798293879538e-07
OrderedDict([('P_impinge_mphys', 1.3868290880508479e-08)])

cdv w.r.t. T:

-2.088757335497182e-07
OrderedDict([('T_impinge_mphys', 5.804483645738395e-12)])

Code versions

• Operating System: Linux 5.4.72-microsoft-standard-WSL2 x86_64
• Python: 3.8
• OpenMPI: 4.1.1
• CGNS: 4.1.2
• PETSc: 3.15.5
• This software: Current Release, but issue occurs with previous releases, notably 2.30 with PETSc: 3.12.5

[sseraj]

The beta derivative is wrong because of a bug in ADflow. The value needs to be converted from 1/radians to 1/degrees as done for alpha:

adflow/adflow/pyADflow.py

Lines 3599 to 3600 in 196caaa

	if key == "alpha":
	funcsSens[dvName] *= numpy.pi / 180.0

Did you use a step size of 1e-5 for all the DVs? If so, I suspect the step size is too small for P and T. I would recommend doing a step size study or comparing against complex step.

[garobed1]

I see, the 180/pi factor does correct the beta derivative. I'll correct for that in my work. Given that, I do suspect that the other derivatives are a case of the step size being too small, but I'll have to check tomorrow. Thanks!

[garobed1]

While the beta derivative is resolved, I'm still getting disagreeing results for pressure and temperature. To isolate the issue from my own work, I tried the example NACA0012 problem with RANS-SA and got similarly disagreeing results between analytic and central-difference derivatives for the pressure variable (off by 1 order of magnitude and sign). Interestingly, if I define the Aeroproblem using altitude, gradients with respect to altitude agree with each other, but if I take the pressure and temperature defined by that same altitude and define the AeroProblem that way, the analytic derivative with respect to pressure seems incorrect.

-4.046138085494165e-10 #central difference h=1e-3
{'fc_cdv': OrderedDict([('P_fc', 7.864265527134046e-11)])} #analytic

I tried complex step (h=1e-30) through ADFlow and got a value that agrees with the central difference

-4.012000967017802e-10 #complex step

Code for this:

import numpy
from mpi4py import MPI
from baseclasses import *
from adflow import ADFLOW, ADFLOW_C
from baseclasses.problems.ICAOAtmosphere import ICAOAtmosphere

# ======================================================================
#         Input Information -- Modify accordingly!
# ======================================================================
outputDirectory = './'
gridFile = './adflow/inputFiles/input_files/naca0012_rans-L2.cgns'
alpha = 2.0
mach = 1.78
areaRef = 45.5
chordRef = 3.25
MGCycle = 'sg'
altitude = 10000
name = 'fc'

aeroOptions = {
# Common Parameters
'gridFile':gridFile,
'outputDirectory':outputDirectory,
'printIterations':True,

# Physics Parameters
'equationType':'RANS',
'turbulenceModel':'SA',

# Common Parameters
'CFL':1.5,
'CFLCoarse':1.25,
'MGCycle':MGCycle,
'MGStartLevel':-1,
'nCyclesCoarse':250,
'nCycles':1000000,
'monitorvariables':['resrho','cl','cd'],
'useNKSolver':True,
'NKSwitchTol':1e-4,

# Convergence Parameters
'L2Convergence':1e-12,
'L2ConvergenceCoarse':1e-2,
}


# Aerodynamic problem description
atm = ICAOAtmosphere()
P, T = atm(altitude)
ap = AeroProblem(name=name, alpha=alpha, mach=mach, P=P, T=T, #altitude=altitude,
areaRef=areaRef, chordRef=chordRef,
evalFuncs=['cdv'])
ap.addDV('P')
# Create solver
#CFDSolver = ADFLOW(options=aeroOptions)
CFDSolver = ADFLOW_C(options=aeroOptions)


# Solve and evaluate functions
# funcs = {}
# funcsens = {}
# CFDSolver(ap)
# CFDSolver.evalFunctions(ap, funcs)
# CFDSolver.evalFunctionsSens(ap, funcsens)

# central difference check, directional derivative
# step size
h = 1e-40

# just do full fd/cs
fd = 0.


funcs2 = {}
ap.P += h*1j
CFDSolver(ap)
CFDSolver.evalFunctions(ap, funcs2)

# funcs2 = {}
# ap.P += h
# CFDSolver(ap)
# CFDSolver.evalFunctions(ap, funcs2)

# funcs3 = {}
# ap.P -= 2*h
# CFDSolver(ap)
# CFDSolver.evalFunctions(ap, funcs3)
# fd = (funcs2['fc_cdv'] - funcs3['fc_cdv'])/(2*h)
cs = numpy.imag(funcs2['fc_cdv'])/h

if MPI.COMM_WORLD.rank == 0:
    print(cs)
    # print(fd)
    # print(funcsens)

[sseraj]

Thanks for the detailed report. I am able to reproduce the error. It looks like a bug, but I don't have a good idea on how to fix it right now. Please update the issue if you find any leads.

[garobed1]

I got some time to investigate this myself a bit, but I haven't been able to find anything conclusive. The only insight I have is that I'm inclined to trust the AD pressure derivative because the altitude derivative depends on it in the chain rule, and that derivative is correct with the value given by the AD pressure derivative.

From pyADFlow:

if key == 'altitude':
               # This design variable is special. It combines changes
               # in temperature, pressure and density into a single
               # variable. Since we have derivatives for T, P and
               # rho, we simply chain rule it back to the the
               # altitude variable.
               self.curAP.evalFunctionsSens(tmp, ['P', 'T', 'rho'])

               # Extract the derivatives wrt the independent
               # parameters in ADflow
               dIdP = dIda[self.possibleAeroDVs['p']]  #AD pressure derivative
               dIdT = dIda[self.possibleAeroDVs['t']]
               dIdrho = dIda[self.possibleAeroDVs['rho']]

               # Chain-rule to get the final derivative:
               funcsSens[dvName] = (
                   tmp[self.curAP['P']][dvName]*dIdP + 
                   tmp[self.curAP['T']][dvName]*dIdT +
                   tmp[self.curAP['rho']][dvName]*dIdrho)

tmp[self.curAP['P']][dvName] = -4.042423542031349 #pressure derivative with respect to altitude
dIdP = 7.34674202076354e-11 #AD pressure derivative

dIdPfd =  -3.8665143135829805e-10 #FD pressure derivative

1.2145280795286122e-09 #AD altitude derivative with AD pressure derivative
3.07452336e-9                   #AD altitude derivative if we used the FD pressure derivative instead
1.2100201483150608e-09 #FD altitude derivative

If the altitude derivative used the finite-differenced/complex-step pressure derivative, it would be way off. This still doesn't explain why the AD pressure derivative disagrees with finite difference, but it lends credibility to the idea that the pressure derivative is actually correct.

I don't have any further insight, however.