Moved the application of Sobolev gradient smoothing to SU2_DOT_AD.

 - Added the necessary solver, numerics, and variable classes there.
 - The method is called as an alternative parallel to SetProject_AD.

Common/include/option_structure.hpp

@@ -533,8 +533,8 @@ const int EL_TETRA = 0;   /*!< \brief Elements of four nodes (3D). */
 const int EL_HEXA  = 1;   /*!< \brief Elements of eight nodes (3D). */
 const int EL_PYRAM = 2;   /*!< \brief Elements of five nodes (3D). */
 const int EL_PRISM = 3;   /*!< \brief Elements of six nodes (3D). */
-const int EL_TETRA2 = 4;		/*!< \brief Elements of four nodes, with second order gauss quadrature (3D). */
-const int EL_PYRAM2 = 5;		/*!< \brief Elements of five nodes, with third order gauss quadrature (3D). */
+const int EL_TETRA2 = 4;  /*!< \brief Elements of four nodes, with second order gauss quadrature (3D). */
+const int EL_PYRAM2 = 5;  /*!< \brief Elements of five nodes, with third order gauss quadrature (3D). */
 
 /*!
  * \brief Types of spatial discretizations

Common/src/CConfig.cpp

@@ -5192,7 +5192,7 @@ void CConfig::SetMarkers(SU2_COMPONENT val_software) {
   Marker_All_Turbomachinery       = new unsigned short[nMarker_All] (); // Store whether the boundary is in needed for Turbomachinery computations.
   Marker_All_TurbomachineryFlag   = new unsigned short[nMarker_All] (); // Store whether the boundary has a flag for Turbomachinery computations.
   Marker_All_MixingPlaneInterface = new unsigned short[nMarker_All] (); // Store whether the boundary has a in the MixingPlane interface.
-  Marker_All_SobolevBC      = new unsigned short[nMarker_All] ();	// Store wether the boundary should apply to the gradient smoothing.
+  Marker_All_SobolevBC      = new unsigned short[nMarker_All] (); // Store wether the boundary should apply to the gradient smoothing.
 
   for (iMarker_All = 0; iMarker_All < nMarker_All; iMarker_All++) {
     Marker_All_TagBound[iMarker_All] = "SEND_RECEIVE";

Common/src/geometry/elements/CTETRA4.cpp

@@ -50,10 +50,10 @@ CTETRA4::CTETRA4() : CElementWithKnownSizes<NGAUSS,NNODE,NDIM>() {
     Eta = GaussCoord[iGauss][1];
     Zeta = GaussCoord[iGauss][2];
 
-    val_Ni = Xi;						  GaussPoint[iGauss].SetNi(val_Ni,0);
-    val_Ni = Eta;						  GaussPoint[iGauss].SetNi(val_Ni,1);
-    val_Ni = 1.0-Xi-Eta-Zeta;	GaussPoint[iGauss].SetNi(val_Ni,2);
-    val_Ni = Zeta;					  GaussPoint[iGauss].SetNi(val_Ni,3);
+    val_Ni = Xi;              GaussPoint[iGauss].SetNi(val_Ni,0);
+    val_Ni = Eta;             GaussPoint[iGauss].SetNi(val_Ni,1);
+    val_Ni = 1.0-Xi-Eta-Zeta; GaussPoint[iGauss].SetNi(val_Ni,2);
+    val_Ni = Zeta;            GaussPoint[iGauss].SetNi(val_Ni,3);
 
     /*--- dN/d xi, dN/d eta, dN/d zeta ---*/
 

SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp

@@ -122,11 +122,6 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver {
    */
   void SecondaryRecording(void);
 
-  /*!
-   * \brief Postprocess the calculated derivatives with gradient smoothing
-   */
-  void DerivativeTreatment(void);
-
   /*!
    * \brief gets Convergence on physical time scale, (deactivated in adjoint case)
    * \return false

SU2_CFD/include/numerics/CNumerics.hpp

@@ -35,7 +35,6 @@
 
 #include "../../../Common/include/CConfig.hpp"
 #include "../../../Common/include/linear_algebra/blas_structure.hpp"
-#include "../../../Common/include/containers/C2DContainer.hpp"
 
 class CElement;
 class CFluidModel;
@@ -833,12 +832,6 @@ class CNumerics {
     Coord_j = val_coord_j;
   }
 
-  /*!
-   * \brief A virtual member
-   * \param[in] val_coord - vector of coordinates.
-   */
-  virtual void SetCoord(su2activematrix& val_coord) {}
-
   /*!
    * \brief Set the velocity of the computational grid.
    * \param[in] val_gridvel_i - Grid velocity of the point i.

SU2_CFD/include/solvers/CDiscAdjSolver.hpp

@@ -36,7 +36,7 @@
  * \ingroup Discrete_Adjoint
  * \author T. Albring
  */
-class CDiscAdjSolver : public CSolver {
+class CDiscAdjSolver final : public CSolver {
 protected:
   static constexpr size_t MAXNDIM = 3;  /*!< \brief Max number of space dimensions, used in some static arrays. */
   static constexpr size_t MAXNVAR = 32; /*!< \brief Max number of variables, for static arrays. */

SU2_CFD/include/solvers/CSolver.hpp

@@ -4164,8 +4164,7 @@ class CSolver {
    * \param[in] config - Definition of the particular problem.
    */
   virtual void ApplyGradientSmoothingVolume(CGeometry *geometry,
-                                            CSolver *solver,
-                                            CNumerics **numerics,
+                                            CNumerics *numerics,
                                             CConfig *config) { }
 
   /*!
@@ -4177,20 +4176,19 @@ class CSolver {
    *
    */
   virtual void ApplyGradientSmoothingSurface(CGeometry *geometry,
-                                               CSolver *solver,
-                                               CNumerics **numerics,
-                                               CConfig *config,
-                                               unsigned long val_marker) { }
+                                             CNumerics *numerics,
+                                             CConfig *config,
+                                             unsigned long val_marker) { }
 
   /*!
    * \brief All steps required for smoothing the whole system on DV level in an iterative way
    */
   virtual void ApplyGradientSmoothingDV(CGeometry *geometry,
-                                        CSolver *solver,
-                                        CNumerics **numerics,
+                                        CNumerics *numerics,
                                         CSurfaceMovement *surface_movement,
                                         CVolumetricMovement *grid_movement,
-                                        CConfig *config) { }
+                                        CConfig *config,
+                                        su2double** Gradient) { }
 
   /*!
    * \brief A virtual member.
@@ -4200,27 +4198,29 @@ class CSolver {
   virtual void RecordTapeAndCalculateOriginalGradient(CGeometry *geometry,
                                                      CSurfaceMovement *surface_movement,
                                                      CVolumetricMovement *grid_movement,
-                                                     CConfig *config) { }
+                                                     CConfig *config,
+                                                     su2double** Gradient) { }
 
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver - The solver container holding all terms of the solution.
    * \param[in] config - Definition of the particular problem.
    */
-  virtual void OutputSensitivity(CGeometry *geometry,
-                                 CConfig *config,
-                                 CSolver *solver) { }
+  virtual void ReadSens2Geometry(CGeometry *geometry,
+                                 CConfig *config) { }
 
   /*!
    * \brief A virtual member.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
    */
-  inline virtual vector<su2double> GetDeltaP() { return vector<su2double>(); }
+  virtual void WriteSens2Geometry(CGeometry *geometry,
+                                 CConfig *config) { }
 
   /*!
    * \brief A virtual member.
    */
-  virtual void GetHessianMatrix() {}
+  inline virtual vector<su2double> GetDeltaP() { return vector<su2double>(); }
 
   /*!
    * \brief Routine that sets the flag controlling implicit treatment for periodic BCs.

SU2_CFD/include/solvers/CSolverFactory.hpp

@@ -61,7 +61,6 @@ enum class SUB_SOLVER_TYPE {
   MESH,                    /*!< \brief Mesh solver */
   RADIATION,               /*!< \brief Radiation solver */
   DISC_ADJ_RADIATION,      /*!< \brief Discrete adjoint radiation solver */
-  GRAD_SMOOTH,             /*!< \brief Gradient Smoothing solver for discrete adjoints. */
   NONE
 };
 

SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp

@@ -510,60 +510,4 @@ void CDiscAdjSinglezoneDriver::SecondaryRecording(){
 
   AD::ClearAdjoints();
 
-  /*--- If necessary smooth the calculated geometry sensitivities ---*/
-
-  if (config->GetSmoothGradient()) {
-    DerivativeTreatment();
-  }
-
-}
-
-/*--- Main routine to call the Sobolev smoothing solver and postprocess the derivatives. ---*/
-
-void CDiscAdjSinglezoneDriver::DerivativeTreatment() {
-
-  if (rank == MASTER_NODE)  cout << "Sobolev Smoothing of derivatives is active." << endl;
-
-  /*--- Get the sensitivities from the adjoint solver to work with. ---*/
-  solver[GRADIENT_SMOOTHING]->SetSensitivity(geometry,config,solver[MainSolver]);
-
-  /*--- Apply the smoothing procedure on the mesh level. ---*/
-  if (config->GetSobMode()==MESH_LEVEL || config->GetSobMode()==DEBUG ) {
-
-    if (rank == MASTER_NODE)  cout << "  working on mesh level (including debug mode)" << endl;
-
-    /*--- Work with the surface derivatives. ---*/
-    if (config->GetSmoothOnSurface()) {
-
-      /*--- Select DV marker, or NOT_AVAILABLE if none is specified. ---*/
-      unsigned long  dvMarker = BC_TYPE::NOT_AVAILABLE;
-      for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
-        if ( config->GetMarker_All_DV(iMarker) == YES ) {
-          dvMarker = iMarker;
-        }
-      }
-      solver[GRADIENT_SMOOTHING]->ApplyGradientSmoothingSurface(geometry, solver[MainSolver], numerics[GRADIENT_SMOOTHING], config, dvMarker);
-
-    /*--- Work with the volume derivatives. ---*/
-    } else {
-      solver[GRADIENT_SMOOTHING]->ApplyGradientSmoothingVolume(geometry, solver[MainSolver], numerics[GRADIENT_SMOOTHING], config);
-    }
-
-    /*--- After appling the solver write the results back ---*/
-    solver[GRADIENT_SMOOTHING]->OutputSensitivity(geometry,config,solver[ADJFLOW_SOL]);
-
-  /*--- Apply the smoothing procedure on the DV level. ---*/
-  } else if (config->GetSobMode()==PARAM_LEVEL_COMPLETE) {
-
-    solver[GRADIENT_SMOOTHING]->ApplyGradientSmoothingDV(geometry, solver[MainSolver], numerics[GRADIENT_SMOOTHING], surface_movement[ZONE_0], grid_movement[ZONE_0][INST_0], config);
-
-  /*--- For some application, e.g. OneShot, we might only need the original gradient. ---*/
-  } else if (config->GetSobMode()==ONLY_GRAD) {
-    solver[GRADIENT_SMOOTHING]->RecordTapeAndCalculateOriginalGradient(geometry, surface_movement[ZONE_0], grid_movement[ZONE_0][INST_0], config);
-
-  /*--- Warning if choose mode is unsupported. ---*/
-  } else if (config->GetSobMode()==NO_MODUS) {
-    if (rank == MASTER_NODE)  cout << "Unsupported operation modus for the Sobolev Smoothing Solver." << endl;
-  }
-
 }

SU2_CFD/src/drivers/CDriver.cpp

@@ -89,7 +89,6 @@
 #include "../../include/numerics/elasticity/CFEALinearElasticity.hpp"
 #include "../../include/numerics/elasticity/CFEANonlinearElasticity.hpp"
 #include "../../include/numerics/elasticity/nonlinear_models.hpp"
-#include "../../include/numerics/CGradSmoothing.hpp"
 
 #include "../../include/integration/CIntegrationFactory.hpp"
 
@@ -232,14 +231,6 @@ CDriver::CDriver(char* confFile, unsigned short val_nZone, SU2_Comm MPICommunica
 
   }
 
-  /*--- For gradient smoothing we need to initialize the mesh deformation classes.
-   * This is in order to calculate the projections between volume, surface and design parameters ---*/
-  if ( config_container[ZONE_0]->GetSmoothGradient() ) {
-    grid_movement[ZONE_0][INST_0] = new CVolumetricMovement(geometry_container[ZONE_0][INST_0][MESH_0], config_container[ZONE_0]);
-    surface_movement[ZONE_0] = new CSurfaceMovement();
-    surface_movement[ZONE_0]->CopyBoundary(geometry_container[ZONE_0][INST_0][MESH_0], config_container[ZONE_0]);
-  }
-
   /*! --- Compute the wall distance again to correctly compute the derivatives if we are running direct diff mode --- */
   if (driver_config->GetDirectDiff() == D_DESIGN){
     CGeometry::ComputeWallDistance(config_container, geometry_container);
@@ -2259,15 +2250,6 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CGeometry **geometry, CSol
 
   } // end "per-thread" allocation loop
 
-  /*--- If we want to apply the gradient smoothing solver we must initialize the numerics classes. ---*/
-  if(config->GetSmoothGradient()) {
-    if (config->GetSmoothOnSurface()) {
-      numerics[MESH_0][GRADIENT_SMOOTHING][GRAD_TERM] = new CGradSmoothing(nDim-1, config);
-    } else {
-      numerics[MESH_0][GRADIENT_SMOOTHING][GRAD_TERM] = new CGradSmoothing(nDim, config);
-    }
-  }
-
 }
 
 void CDriver::Numerics_Postprocessing(CNumerics *****numerics, CSolver***, CGeometry**,

SU2_CFD/src/meson.build

@@ -77,8 +77,7 @@ su2_cfd_src += files(['variables/CIncNSVariable.cpp',
                       'variables/CIncEulerVariable.cpp',
                       'variables/CEulerVariable.cpp',
                       'variables/CNEMOEulerVariable.cpp',
-                      'variables/CNEMONSVariable.cpp',
-                      'variables/CSobolevSmoothingVariable.cpp'])
+                      'variables/CNEMONSVariable.cpp'])
 
 su2_cfd_src += files(['solvers/CSolverFactory.cpp',
                       'solvers/CAdjEulerSolver.cpp',
@@ -107,8 +106,7 @@ su2_cfd_src += files(['solvers/CSolverFactory.cpp',
                       'solvers/CTransLMSolver.cpp',
                       'solvers/CTurbSolver.cpp',
                       'solvers/CTurbSASolver.cpp',
-                      'solvers/CTurbSSTSolver.cpp',
-                      'solvers/CGradientSmoothingSolver.cpp'])
+                      'solvers/CTurbSSTSolver.cpp'])
 
 su2_cfd_src += files(['numerics/CNumerics.cpp',
                       'numerics/template.cpp',
@@ -141,8 +139,7 @@ su2_cfd_src += files(['numerics/CNumerics.cpp',
                       'numerics/elasticity/CFEAElasticity.cpp',
                       'numerics/elasticity/CFEALinearElasticity.cpp',
                       'numerics/elasticity/CFEANonlinearElasticity.cpp',
-                      'numerics/elasticity/nonlinear_models.cpp',
-                      'numerics/CGradSmoothing.cpp'])
+                      'numerics/elasticity/nonlinear_models.cpp'])
 
 su2_cfd_src += files(['../include/numerics_simd/CNumericsSIMD.cpp'])
 

SU2_CFD/src/solvers/CSolverFactory.cpp

@@ -51,7 +51,6 @@
 #include "../../include/solvers/CBaselineSolver.hpp"
 #include "../../include/solvers/CBaselineSolver_FEM.hpp"
 #include "../../include/solvers/CRadP1Solver.hpp"
-#include "../../include/solvers/CGradientSmoothingSolver.hpp"
 
 map<const CSolver*, SolverMetaData> CSolverFactory::allocatedSolvers;
 
@@ -116,26 +115,22 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon
     case MAIN_SOLVER::DISC_ADJ_EULER:
       solver[FLOW_SOL]    = CreateSubSolver(SUB_SOLVER_TYPE::EULER, solver, geometry, config, iMGLevel);
       solver[ADJFLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_FLOW, solver, geometry, config, iMGLevel);
-      solver[GRADIENT_SMOOTHING] = CreateSubSolver(SUB_SOLVER_TYPE::GRAD_SMOOTH, solver, geometry, config, iMGLevel);
       break;
     case MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES:
       solver[FLOW_SOL]    = CreateSubSolver(SUB_SOLVER_TYPE::NAVIER_STOKES, solver, geometry, config, iMGLevel);
       solver[ADJFLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_FLOW, solver, geometry, config, iMGLevel);
-      solver[GRADIENT_SMOOTHING] = CreateSubSolver(SUB_SOLVER_TYPE::GRAD_SMOOTH, solver, geometry, config, iMGLevel);
       break;
     case MAIN_SOLVER::DISC_ADJ_RANS:
       solver[FLOW_SOL]    = CreateSubSolver(SUB_SOLVER_TYPE::NAVIER_STOKES, solver, geometry, config, iMGLevel);
       solver[ADJFLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_FLOW, solver, geometry, config, iMGLevel);
       solver[TURB_SOL]    = CreateSubSolver(SUB_SOLVER_TYPE::TURB, solver, geometry, config, iMGLevel);
       solver[ADJTURB_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_TURB, solver, geometry, config, iMGLevel);
-      solver[GRADIENT_SMOOTHING] = CreateSubSolver(SUB_SOLVER_TYPE::GRAD_SMOOTH, solver, geometry, config, iMGLevel);
       break;
     case MAIN_SOLVER::DISC_ADJ_INC_EULER:
       solver[FLOW_SOL]    = CreateSubSolver(SUB_SOLVER_TYPE::INC_EULER, solver, geometry, config, iMGLevel);
       solver[ADJFLOW_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_FLOW, solver, geometry, config, iMGLevel);
       solver[RAD_SOL]     = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel);
       solver[ADJRAD_SOL]  = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_RADIATION, solver, geometry, config, iMGLevel);
-      solver[GRADIENT_SMOOTHING] = CreateSubSolver(SUB_SOLVER_TYPE::GRAD_SMOOTH, solver, geometry, config, iMGLevel);
       break;
     case MAIN_SOLVER::DISC_ADJ_INC_NAVIER_STOKES:
       solver[FLOW_SOL]    = CreateSubSolver(SUB_SOLVER_TYPE::INC_NAVIER_STOKES, solver, geometry, config, iMGLevel);
@@ -144,7 +139,6 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon
       solver[ADJHEAT_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_HEAT, solver, geometry, config, iMGLevel);
       solver[RAD_SOL]     = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel);
       solver[ADJRAD_SOL]  = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_RADIATION, solver, geometry, config, iMGLevel);
-      solver[GRADIENT_SMOOTHING] = CreateSubSolver(SUB_SOLVER_TYPE::GRAD_SMOOTH, solver, geometry, config, iMGLevel);
       break;
     case MAIN_SOLVER::DISC_ADJ_INC_RANS:
       solver[FLOW_SOL]    = CreateSubSolver(SUB_SOLVER_TYPE::INC_NAVIER_STOKES, solver, geometry, config, iMGLevel);
@@ -155,7 +149,6 @@ CSolver** CSolverFactory::CreateSolverContainer(MAIN_SOLVER kindMainSolver, CCon
       solver[ADJTURB_SOL] = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_TURB, solver, geometry, config, iMGLevel);
       solver[RAD_SOL]     = CreateSubSolver(SUB_SOLVER_TYPE::RADIATION, solver, geometry, config, iMGLevel);
       solver[ADJRAD_SOL]  = CreateSubSolver(SUB_SOLVER_TYPE::DISC_ADJ_RADIATION, solver, geometry, config, iMGLevel);
-      solver[GRADIENT_SMOOTHING] = CreateSubSolver(SUB_SOLVER_TYPE::GRAD_SMOOTH, solver, geometry, config, iMGLevel);
     break;
     case MAIN_SOLVER::DISC_ADJ_HEAT:
       solver[HEAT_SOL]    = CreateSubSolver(SUB_SOLVER_TYPE::HEAT, solver, geometry, config, iMGLevel);
@@ -305,11 +298,6 @@ CSolver* CSolverFactory::CreateSubSolver(SUB_SOLVER_TYPE kindSolver, CSolver **s
       }
       metaData.integrationType = INTEGRATION_TYPE::DEFAULT;
       break;
-    case SUB_SOLVER_TYPE::GRAD_SMOOTH:
-      if (iMGLevel == MESH_0 && config->GetSmoothGradient() ) {
-        genericSolver = new CGradientSmoothingSolver(geometry, config);
-      }
-      break;
     default:
       SU2_MPI::Error("No proper allocation found for requested sub solver", CURRENT_FUNCTION);
       break;

SU2_DOT/include/SU2_DOT.hpp

@@ -46,6 +46,9 @@
 #include "../../SU2_CFD/include/output/CBaselineOutput.hpp"
 #include "../../SU2_CFD/include/solvers/CBaselineSolver.hpp"
 
+#include "solvers/CGradientSmoothingSolver.hpp"
+#include "numerics/CGradSmoothing.hpp"
+
 using namespace std;
 
 
@@ -87,3 +90,13 @@ void OutputGradient(su2double** Gradient, CConfig* config, ofstream& Gradient_fi
  */
 
 void SetSensitivity_Files(CGeometry ***geometry, CConfig **config, unsigned short val_nZone);
+
+/*!
+ * \brief Treatment of derivatives with the Sobolev smoothing solver.
+ * \param[in] geometry - Geometrical definition of the problem.
+ * \param[in] config - Definition of the particular problem.
+ * \param[in] surface_movement - Surface movement class of the problem.
+ * \param[in] Gradient_file - Output file to store the gradient data.
+ */
+
+void DerivativeTreatment(CGeometry *geometry, CConfig *config, CVolumetricMovement *grid_movement, CSurfaceMovement *surface_movement, su2double **Gradient);

SU2_CFD/include/numerics/CGradSmoothing.hpp → SU2_DOT/include/numerics/CGradSmoothing.hpp

@@ -27,8 +27,8 @@
 
 #pragma once
 
-#include "CNumerics.hpp"
 #include "../../../Common/include/geometry/elements/CElement.hpp"
+#include "../../../SU2_CFD/include/numerics/CNumerics.hpp"
 
 /*!
  * \class CGradSmoothing