Fix preaccumulation error, improve recording of RANS solvers

SU2_CFD/include/solvers/CTurbSolver.hpp

@@ -140,4 +140,11 @@ class CTurbSolver : public CScalarSolver<CTurbVariable> {
       Inlet_TurbVars[val_marker][val_vertex][val_dim] = val_turb_var;
   }
 
+  /*!
+   * \brief Register additional In- or Outputs for RANS.
+   * \param[in] input - Boolean whether In- or Output should be registered.
+   * \param[in] config - The particular config.
+   * \returns The number of extra variables.
+   */
+  unsigned long RegisterSolutionExtra(bool input, const CConfig* config) final;
 };

SU2_CFD/include/variables/CTurbVariable.hpp

@@ -100,5 +100,9 @@ class CTurbVariable : public CScalarVariable {
    */
   inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; }
 
+  /*!
+   * \brief Register eddy viscosity (muT) as Input or Output of an AD recording.
+   * \param[in] input - Boolean whether In- or Output should be registered.
+   */
+  void RegisterEddyViscosity(bool input);
 };
-

SU2_CFD/include/variables/CVariable.hpp

@@ -137,6 +137,20 @@ class CVariable {
     RegisterContainer(input, variable, &ad_index);
   }
 
+  void RegisterContainer(bool input, su2activevector& variable, su2vector<int>* ad_index = nullptr) {
+    const auto nPoint = variable.rows();
+    SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads()))
+    for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
+
+      if (input) AD::RegisterInput(variable(iPoint));
+      else AD::RegisterOutput(variable(iPoint));
+
+      if (ad_index) AD::SetIndex((*ad_index)(iPoint), variable(iPoint));
+
+    }
+    END_SU2_OMP_FOR
+  }
+
 public:
   /*--- Disable copy and assignment. ---*/
   CVariable(const CVariable&) = delete;

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

@@ -162,6 +162,7 @@ void CDiscAdjSolver::RegisterSolution(CGeometry *geometry, CConfig *config) {
   /*--- Boolean true indicates that an input is registered ---*/
   direct_solver->GetNodes()->RegisterSolution(true);
 
+  /*--- Register quantities that are no solver variables but further inputs/outputs of the (outer) iteration. ---*/
   direct_solver->RegisterSolutionExtra(true, config);
 
   if (time_n_needed)

SU2_CFD/src/solvers/CTurbSolver.cpp

@@ -236,3 +236,12 @@ void CTurbSolver::Impose_Fixed_Values(const CGeometry *geometry, const CConfig *
   }
 
 }
+
+unsigned long CTurbSolver::RegisterSolutionExtra(bool input, const CConfig* config) {
+
+  /*--- Register muT as input/output of a RANS iteration. ---*/
+  nodes->RegisterEddyViscosity(input);
+
+  /*--- We don't need to save adjoint values for muT. ---*/
+  return 0;
+}

SU2_CFD/src/variables/CNSVariable.cpp

@@ -114,6 +114,8 @@ void CNSVariable::SetRoe_Dissipation_FD(unsigned long iPoint, su2double val_wall
   AD::SetPreaccIn(Primitive(iPoint, indices.EddyViscosity()));
   /*--- Laminar viscosity --- */
   AD::SetPreaccIn(Primitive(iPoint, indices.LaminarViscosity()));
+  /*--- Density; GetDensity reads from Solution (not Primitive) at index 0 ---*/
+  AD::SetPreaccIn(Solution(iPoint, 0));
 
   su2double uijuij = 0.0;
 

SU2_CFD/src/variables/CTurbVariable.cpp

@@ -36,3 +36,7 @@ CTurbVariable::CTurbVariable(unsigned long npoint, unsigned long ndim, unsigned
     intermittency.resize(nPoint) = su2double(1.0);
 
    }
+
+void CTurbVariable::RegisterEddyViscosity(bool input) {
+  RegisterContainer(input, muT);
+}