From aee0c30a3edb9856ce2fc6f7a92772ba9a80aae2 Mon Sep 17 00:00:00 2001
From: Vaggelis Papoutsis <vaggelisp@gmail.com>
Date: Tue, 1 Dec 2020 17:34:32 +0200
Subject: [PATCH] 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.
---
 .../sensitivitySurfaceIncompressible.C        | 20 +++++++++++--------
 .../sensitivitySurfaceIncompressible.H        |  7 +++++--
 .../sensitivitySurfacePointsIncompressible.C  | 19 ++++++++++--------
 .../sensitivitySurfacePointsIncompressible.H  |  3 +++
 4 files changed, 31 insertions(+), 18 deletions(-)

diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C
index 0f111458949..69aa297c62f 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.C
+++ b/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_)
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H
index 59cc5752916..bbec0abfa1d 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurface/sensitivitySurfaceIncompressible.H
+++ b/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_;
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C
index ce9654cf16a..484c3f889c2 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.C
+++ b/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_)
diff --git a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.H b/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.H
index bf889a0f4d5..abe1a2e5456 100644
--- a/src/optimisation/adjointOptimisation/adjoint/optimisation/adjointSensitivity/incompressible/sensitivitySurfacePoints/sensitivitySurfacePointsIncompressible.H
+++ b/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_;
 
-- 
GitLab