Incorrect AP variable function sensitivities for simple 2D shock problem

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

Start with the provided Python reproducer and the ADFLOW evalFunctionsSens entry point, then compare its sensitivities against the central-difference checks for beta, P, T, and mach. The CGNS mesh is not attached, so obtain a reproducible mesh or isolate the setup that triggers the discrepancy. Done means explaining and correcting the inconsistent sensitivities, with the mach result still matching the finite-difference check.