Improve the Euler wall Jacobian

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1141,59 +1141,26 @@
       for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim] / Area;
     }
 
-    su2double* V_reflected = GetCharacPrimVar(val_marker, iVertex);
-
-    /*--- Grid movement ---*/
-    if (dynamic_grid)
-      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);
-
-    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)
-      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];
-
-    /*--- 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));
-
-    /*--- 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];
-    }
-
     /*--- Explicitly set the velocity components normal to the symmetry plane to zero.
      * This is necessary because the modification of the residual leaves the problem
      * underconstrained (the normal residual is zero regardless of the normal velocity). ---*/
 
     su2double* solutionOld = nodes->GetSolution_Old(iPoint);
 
-    su2double gridVel[MAXNVAR] = {};
+    su2double gridVel[MAXNDIM] = {}, qn = 0.0;
     if (dynamic_grid) {
       for (auto iDim = 0u; iDim < nDim; iDim++) {
         gridVel[iDim] = geometry->nodes->GetGridVel(iPoint)[iDim];
+        qn += gridVel[iDim] * Normal[iDim];
       }
       if (FlowRegime == ENUM_REGIME::COMPRESSIBLE) {
         for(auto iDim = 0u; iDim < nDim; iDim++) {
           /*--- Multiply by density since we are correcting conservative variables. ---*/
           gridVel[iDim] *= nodes->GetDensity(iPoint);
         }
+        /*--- Work of pressure forces. This is not needed for incompressible regimes
+         * because we use Cv == Cp. ---*/
+        LinSysRes(iPoint, iVel + nDim) -= qn * nodes->GetPressure(iPoint);
       }
     }
     su2double vp = 0.0;
@@ -1215,7 +1182,30 @@
 
     /*--- 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 (unsigned short iVar = 0; iVar < nVar; iVar++)
+        mat[iVar * nVar + iVar] = 1;
+      for (unsigned short iDim = 0; iDim < nDim; iDim++)
+        for (unsigned short jDim = 0; 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 (unsigned short iVar = 0; iVar < nVar * nVar; iVar++) jac[iVar] = block[iVar];
+
+        CBlasStructure().gemm(nVar, nVar, nVar, mat, jac, newJac, config);
+
+        for (unsigned short iVar = 0; 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. ---*/

TestCases/rotating/naca0012/rot_NACA0012.cfg

@@ -52,7 +52,7 @@ MARKER_DESIGNING = ( airfoil )
 % ------------- COMMON PARAMETERS TO DEFINE THE NUMERICAL METHOD --------------%
 %
 NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
-CFL_NUMBER= 1e4
+CFL_NUMBER= 100
 CFL_ADAPT= NO
 CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
 RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )