ENH: change in the discretization of part of the (E)SI sensitivity terms

Part of the (E)SI shape sensitivities depends of grad(Ua) & nf computed
on the boundary. Up to now, the code was only computing the normal part
of grad(Ua), to avoid the potentially spurious tangential component
which is computed on the cell center and extrapolated to the boundary
faces. However, for some objectives that are strongly related to the
stresses (e.g. moment, stresses), including also the tangential part of
grad(Ua) is necessary for E-SI to replicate the outcome of FI.

Extensive testing on a number of objectives/cases showed
- No regression when including the tangential part
- Improved behaviour in some rare cases (moment, stresses)

Hence, the tangential part is now included by default. The previous code
behaviour can be replicated by setting the useSnGradInTranposeStresses
flag to true.

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C

@@ -240,6 +240,7 @@ sensitivitySurface::sensitivitySurface
     includePressureTerm_(false),
     includeGradStressTerm_(false),
     includeTransposeStresses_(false),
+    useSnGradInTranposeStresses_(false),
     includeDivTerm_(false),
     includeDistance_(false),
     includeMeshMovement_(false),
@@ -332,6 +333,8 @@ void sensitivitySurface::read()
         dict().getOrDefault<bool>("includeGradStressTerm", true);
     includeTransposeStresses_ =
         dict().getOrDefault<bool>("includeTransposeStresses", true);
+    useSnGradInTranposeStresses_ =
+        dict().getOrDefault<bool>("useSnGradInTranposeStresses", false);
     includeDivTerm_ = dict().getOrDefault<bool>("includeDivTerm", false);
     includeDistance_ =
         dict().getOrDefault<bool>
@@ -528,16 +531,17 @@ void sensitivitySurface::accumulateIntegrand(const scalar dt)
 
         if (includeTransposeStresses_)
         {
+            vectorField gradUaNf
+                (
+                    useSnGradInTranposeStresses_ ?
+                    (Ua.boundaryField()[patchI].snGrad() & nf)*nf :
+                    (gradUa.boundaryField()[patchI] & nf)
+                );
+
             stressTerm -=
                 nuEff.boundaryField()[patchI]
-              * (
-                    // Note: in case of laminar or low-Re flows,
-                    // includes a spurious tangential gradUa component
-                    // (gradUa.boundaryField()[patchI] & nf)
-                    ((Ua.boundaryField()[patchI].snGrad() &nf)*nf)
-                    & U.boundaryField()[patchI].snGrad()
-                )
-              * nf;
+               *(gradUaNf & U.boundaryField()[patchI].snGrad())
+               *nf;
         }
 
         if (includeDivTerm_)

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2007-2019 PCOpt/NTUA
-    Copyright (C) 2013-2019 FOSS GP
+    Copyright (C) 2007-2020 PCOpt/NTUA
+    Copyright (C) 2013-2020 FOSS GP
     Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -79,6 +79,9 @@ protected:
         //- Include the transpose part of the adjoint stresses
         bool includeTransposeStresses_;
 
+        //- Use snGrad in the transpose part of the adjoint stresses
+        bool useSnGradInTranposeStresses_;
+
         //- Include the term from the deviatoric part of the stresses
         bool includeDivTerm_;
 

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C

@@ -61,6 +61,8 @@ void sensitivitySurfacePoints::read()
         dict().getOrDefault<bool>("includeGradStressTerm", true);
     includeTransposeStresses_ =
         dict().getOrDefault<bool>("includeTransposeStresses", true);
+    useSnGradInTranposeStresses_ =
+        dict().getOrDefault<bool>("useSnGradInTranposeStresses", false);
     includeDivTerm_ =
         dict().getOrDefault<bool>("includeDivTerm", false);
     includeDistance_ =
@@ -351,6 +353,7 @@ sensitivitySurfacePoints::sensitivitySurfacePoints
     includePressureTerm_(false),
     includeGradStressTerm_(false),
     includeTransposeStresses_(false),
+    useSnGradInTranposeStresses_(false),
     includeDivTerm_(false),
     includeDistance_(false),
     includeMeshMovement_(false),
@@ -544,16 +547,16 @@ void sensitivitySurfacePoints::accumulateIntegrand(const scalar dt)
 
         if (includeTransposeStresses_)
         {
+            vectorField gradUaNf
+                (
+                    useSnGradInTranposeStresses_ ?
+                    (Ua.boundaryField()[patchI].snGrad() & nf)*nf :
+                    (gradUa.boundaryField()[patchI] & nf)
+                );
             stressTerm -=
                 nuEff.boundaryField()[patchI]
-               *(
-                    // Note: in case of laminar or low-Re flows,
-                    // includes a spurious tangential gradUa component
-                    // (gradUa.boundaryField()[patchI] & nf)
-                    ((Ua.boundaryField()[patchI].snGrad() &nf)*nf)
-                  & U.boundaryField()[patchI].snGrad()
-                )
-              * nf;
+               *(gradUaNf & U.boundaryField()[patchI].snGrad())
+               *nf;
         }
 
         if (includeDivTerm_)

src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.H

@@ -79,6 +79,9 @@ protected:
         //- Include the transpose part of the adjoint stresses
         bool includeTransposeStresses_;
 
+        //- Use snGrad in the transpose part of the adjoint stresses
+        bool useSnGradInTranposeStresses_;
+
         //- Include the term from the deviatoric part of the stresses
         bool includeDivTerm_;