Merge branch 'develop' into feature_custom_source

Author: bigfooted

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1143,41 +1143,54 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
       for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim] / Area;
     }
 
-    su2double* V_reflected = GetCharacPrimVar(val_marker, iVertex);
+    /*--- Energy terms due to grid movement (aka work of pressure forces). ---*/
+    if (dynamic_grid) {
+      su2double* V_reflected = GetCharacPrimVar(val_marker, iVertex);
 
-    /*--- Grid movement ---*/
-    if (dynamic_grid)
-      conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint));
+      conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint),
+                                geometry->nodes->GetGridVel(iPoint));
 
-    /*--- Normal vector for this vertex (negate for outward convention). ---*/
-    for (auto iDim = 0u; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim];
-    conv_numerics->SetNormal(Normal);
+      /*--- Normal vector for this vertex (negate for outward convention). ---*/
+      for (auto iDim = 0u; iDim < nDim; iDim++) Normal[iDim] = -Normal[iDim];
+      conv_numerics->SetNormal(Normal);
 
-    for (auto iVar = 0u; iVar < nPrimVar; iVar++)
-      V_reflected[iVar] = nodes->GetPrimitive(iPoint, iVar);
+      for (auto iVar = 0u; iVar < nPrimVar; iVar++)
+        V_reflected[iVar] = nodes->GetPrimitive(iPoint, iVar);
 
-    su2double ProjVelocity_i = nodes->GetProjVel(iPoint, UnitNormal);
-    /*--- Adjustment to v.n due to grid movement. ---*/
-    if (dynamic_grid)
+      su2double ProjVelocity_i = nodes->GetProjVel(iPoint, UnitNormal);
+      /*--- Adjustment to v.n due to grid movement. ---*/
       ProjVelocity_i -= GeometryToolbox::DotProduct(nDim, geometry->nodes->GetGridVel(iPoint), UnitNormal);
 
-    for (auto iDim = 0u; iDim < nDim; iDim++)
-      V_reflected[iDim + iVel] = nodes->GetVelocity(iPoint, iDim) - ProjVelocity_i * UnitNormal[iDim];
+      for (auto iDim = 0u; iDim < nDim; iDim++)
+        V_reflected[iDim + iVel] = nodes->GetVelocity(iPoint, iDim) - ProjVelocity_i * UnitNormal[iDim];
 
-    /*--- Get current solution at this boundary node ---*/
-    const su2double* V_domain = nodes->GetPrimitive(iPoint);
+      /*--- Get current solution at this boundary node. ---*/
+      const su2double* V_domain = nodes->GetPrimitive(iPoint);
 
-    /*--- Set Primitive and Secondary for numerics class. ---*/
-    conv_numerics->SetPrimitive(V_domain, V_reflected);
-    conv_numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint));
+      /*--- Set Primitive and Secondary for numerics class. ---*/
+      conv_numerics->SetPrimitive(V_domain, V_reflected);
+      conv_numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint));
 
-    /*--- Compute the residual using an upwind scheme. ---*/
-    auto residual = conv_numerics->ComputeResidual(config);
+      /*--- Compute the residual using an upwind scheme. ---*/
+      auto residual = conv_numerics->ComputeResidual(config);
 
-    /*--- We include an update of the continuity and energy here, this is important for stability since
-     * these fluxes include numerical diffusion. ---*/
-    for (auto iVar = 0u; iVar < nVar; iVar++) {
-      if (iVar < iVel || iVar >= iVel + nDim) LinSysRes(iPoint, iVar) += residual.residual[iVar];
+      /*--- Use just the energy fluxes to update the residual, adding the others would
+       * increase numerical diffusion which we wish to avoid if possible. ---*/
+      for (auto iVar = iVel + nDim; iVar < nVar; iVar++) {
+        LinSysRes(iPoint, iVar) += residual.residual[iVar];
+      }
+      if (implicit) {
+        auto* block = Jacobian.GetBlock(iPoint, iPoint);
+        /*--- But in the Jacobian we also include the mass flux, this allows some cases with
+         * motion to use larger CFL, for example pywrapper_translating_naca0012. ---*/
+        for (auto iVar = 0u; iVar < nVar; iVar++) {
+          if (iVar < iVel || iVar >= iVel + nDim) {
+            for (auto jVar = 0u; jVar < nVar; jVar++) {
+              block[iVar * nVar + jVar] += SU2_TYPE::GetValue(residual.jacobian_i[iVar][jVar]);
+            }
+          }
+        }
+      }
     }
 
     /*--- Explicitly set the velocity components normal to the symmetry plane to zero.
@@ -1186,7 +1199,7 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
 
     su2double* solutionOld = nodes->GetSolution_Old(iPoint);
 
-    su2double gridVel[MAXNVAR] = {};
+    su2double gridVel[MAXNDIM] = {};
     if (dynamic_grid) {
       for (auto iDim = 0u; iDim < nDim; iDim++) {
         gridVel[iDim] = geometry->nodes->GetGridVel(iPoint)[iDim];
@@ -1217,7 +1230,30 @@ void CFVMFlowSolverBase<V, FlowRegime>::BC_Sym_Plane(CGeometry* geometry, CSolve
 
     /*--- Jacobian contribution for implicit integration. ---*/
     if (implicit) {
-      Jacobian.AddBlock2Diag(iPoint, residual.jacobian_i);
+      /*--- Modify the Jacobians according to the modification of the residual
+       * J_new = (I - n * n^T) * J where n = {0, nx, ny, nz, 0, ...} ---*/
+      su2double mat[MAXNVAR * MAXNVAR] = {};
+
+      for (auto iVar = 0u; iVar < nVar; iVar++)
+        mat[iVar * nVar + iVar] = 1;
+      for (auto iDim = 0u; iDim < nDim; iDim++)
+        for (auto jDim = 0u; jDim < nDim; jDim++)
+          mat[(iDim + iVel) * nVar + jDim + iVel] -= UnitNormal[iDim] * UnitNormal[jDim];
+
+      auto ModifyJacobian = [&](const unsigned long jPoint) {
+        su2double jac[MAXNVAR * MAXNVAR], newJac[MAXNVAR * MAXNVAR];
+        auto* block = Jacobian.GetBlock(iPoint, jPoint);
+        for (auto iVar = 0u; iVar < nVar * nVar; iVar++) jac[iVar] = block[iVar];
+
+        CBlasStructure().gemm(nVar, nVar, nVar, mat, jac, newJac, config);
+
+        for (auto iVar = 0u; iVar < nVar * nVar; iVar++)
+          block[iVar] = SU2_TYPE::GetValue(newJac[iVar]);
+      };
+      ModifyJacobian(iPoint);
+      for (size_t iNeigh = 0; iNeigh < geometry->nodes->GetnPoint(iPoint); ++iNeigh) {
+        ModifyJacobian(geometry->nodes->GetPoint(iPoint, iNeigh));
+      }
     }
 
     /*--- Correction for multigrid. ---*/

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/CSpeciesFlameletSolver.cpp

@@ -77,12 +77,12 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
   auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
 
   /*--- Retrieve spark ignition parameters for spark-type ignition. ---*/
-  if ((flamelet_config_options.ignition_method == FLAMELET_INIT_TYPE::SPARK) && !config->GetRestart()) {
+  if ((flamelet_config_options.ignition_method == FLAMELET_INIT_TYPE::SPARK)) {
     auto spark_init = flamelet_config_options.spark_init;
     spark_iter_start = ceil(spark_init[4]);
     spark_duration = ceil(spark_init[5]);
     unsigned long iter = config->GetMultizone_Problem() ? config->GetOuterIter() : config->GetInnerIter();
-    ignition = ((iter >= spark_iter_start) && (iter <= (spark_iter_start + spark_duration)) && !config->GetRestart());
+    ignition = ((iter >= spark_iter_start) && (iter <= (spark_iter_start + spark_duration)));
   }
 
   SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem);)
@@ -100,7 +100,7 @@ void CSpeciesFlameletSolver::Preprocessing(CGeometry* geometry, CSolver** solver
       /*--- Apply source terms within spark radius. ---*/
       su2double dist_from_center = 0,
                 spark_radius = flamelet_config_options.spark_init[3];
-      dist_from_center = GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(i_point), flamelet_config_options.flame_init.data());
+      dist_from_center = GeometryToolbox::SquaredDistance(nDim, geometry->nodes->GetCoord(i_point), flamelet_config_options.spark_init.data());
       if (dist_from_center < pow(spark_radius,2)) {
         for (auto iVar = 0u; iVar < nVar; iVar++)
           nodes->SetScalarSource(i_point, iVar, nodes->GetScalarSources(i_point)[iVar] + flamelet_config_options.spark_reaction_rates[iVar]);

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);
+}

TestCases/Tutorials/design/Inc_Turbulent_Bend_Wallfunctions/optim.csv.ref

@@ -1,2 +1,2 @@
 ITER,  avg_dp
-1,  63.82289754521027
+1,  64.90301114659223