Improve the Euler wall Jacobian

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -1141,41 +1141,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.
@@ -1184,7 +1197,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];
@@ -1215,7 +1228,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. ---*/

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

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

TestCases/cont_adj_euler/naca0012/of_grad_cd_disc.dat.ref

@@ -1,39 +1,39 @@
 VARIABLES="VARIABLE"      , "GRADIENT"      , "FINDIFF_STEP"  
-          0               , -2656.0         , 0.001           
-          1               , -7290.69        , 0.001           
-          2               , -11252.7        , 0.001           
-          3               , -12531.4        , 0.001           
-          4               , -10923.7        , 0.001           
-          5               , -6826.86        , 0.001           
-          6               , -843.71         , 0.001           
-          7               , 6462.74         , 0.001           
-          8               , 14656.3         , 0.001           
-          9               , 23308.0         , 0.001           
-          10              , 31688.2         , 0.001           
-          11              , 38407.2         , 0.001           
-          12              , 41274.9         , 0.001           
-          13              , 37812.2         , 0.001           
-          14              , 27093.1         , 0.001           
-          15              , 14049.9         , 0.001           
-          16              , 17583.8         , 0.001           
-          17              , 78620.8         , 0.001           
-          18              , 140614.0        , 0.001           
-          19              , -18263.2        , 0.001           
-          20              , -20877.9        , 0.001           
-          21              , -22859.6        , 0.001           
-          22              , -24751.5        , 0.001           
-          23              , -28392.5        , 0.001           
-          24              , -35514.1        , 0.001           
-          25              , -46855.4        , 0.001           
-          26              , -61697.7        , 0.001           
-          27              , -77631.1        , 0.001           
-          28              , -90537.8        , 0.001           
-          29              , -95036.8        , 0.001           
-          30              , -85899.0        , 0.001           
-          31              , -60885.3        , 0.001           
-          32              , -24333.4        , 0.001           
-          33              , 12334.3         , 0.001           
-          34              , 39948.1         , 0.001           
-          35              , 65452.2         , 0.001           
-          36              , 82215.8         , 0.001           
-          37              , -44066.2        , 0.001           
+          0               , 292.459         , 0.001           
+          1               , -8318.5         , 0.001           
+          2               , -16158.4        , 0.001           
+          3               , -21277.2        , 0.001           
+          4               , -23366.9        , 0.001           
+          5               , -22727.6        , 0.001           
+          6               , -19872.3        , 0.001           
+          7               , -15278.1        , 0.001           
+          8               , -9283.99        , 0.001           
+          9               , -2179.87        , 0.001           
+          10              , 5520.63         , 0.001           
+          11              , 12726.0         , 0.001           
+          12              , 17538.0         , 0.001           
+          13              , 17441.1         , 0.001           
+          14              , 10335.6         , 0.001           
+          15              , -2686.37        , 0.001           
+          16              , -10522.5        , 0.001           
+          17              , 24712.2         , 0.001           
+          18              , 166438.0        , 0.001           
+          19              , -15618.6        , 0.001           
+          20              , -14178.7        , 0.001           
+          21              , -12765.5        , 0.001           
+          22              , -12007.2        , 0.001           
+          23              , -13597.5        , 0.001           
+          24              , -19002.9        , 0.001           
+          25              , -28729.6        , 0.001           
+          26              , -41946.6        , 0.001           
+          27              , -56289.3        , 0.001           
+          28              , -67832.1        , 0.001           
+          29              , -71484.0        , 0.001           
+          30              , -62334.5        , 0.001           
+          31              , -38478.2        , 0.001           
+          32              , -4757.34        , 0.001           
+          33              , 26448.5         , 0.001           
+          34              , 45049.5         , 0.001           
+          35              , 60960.9         , 0.001           
+          36              , 83515.9         , 0.001           
+          37              , 8837.4          , 0.001           

TestCases/cont_adj_euler/naca0012/of_grad_directdiff.dat.ref

@@ -1,4 +1,4 @@
 VARIABLES="VARIABLE"      , "DRAG"          , "EFFICIENCY"    , "FORCE_X"       , "FORCE_Y"       , "FORCE_Z"       , "LIFT"          , "MOMENT_X"      , "MOMENT_Y"      , "MOMENT_Z"      , "SIDEFORCE"     
-          0               , 0.2405003164    , -108.3397202    , 0.2752804375    , -1.591326829    , 0.0             , -1.596953347    , 0.0             , 0.0             , 1.18511582      , 0.0             
-          1               , 0.4163088272    , -176.5661566    , 0.4692829784    , -2.423229574    , 0.0             , -2.432890265    , 0.0             , 0.0             , 1.020271816     , 0.0             
-          2               , 0.5558397416    , -231.7769399    , 0.6239540061    , -3.115569035    , 0.0             , -3.128439099    , 0.0             , 0.0             , 0.6093229584    , 0.0             
+          0               , 0.24008251      , -117.3444057    , 0.2742430499    , -1.56293638     , 0.0             , -1.568547024    , 0.0             , 0.0             , 1.189018284     , 0.0             
+          1               , 0.4064005433    , -189.77779      , 0.4586830911    , -2.391641926    , 0.0             , -2.401078899    , 0.0             , 0.0             , 1.030793484     , 0.0             
+          2               , 0.5421052294    , -249.5397676    , 0.6095878335    , -3.086770277    , 0.0             , -3.099333798    , 0.0             , 0.0             , 0.6218682473    , 0.0             

TestCases/cont_adj_euler/wedge/of_grad_combo.dat.ref

@@ -1,5 +1,5 @@
 VARIABLES="VARIABLE"      , "GRADIENT"      , "FINDIFF_STEP"  
-          0               , 0.00770904      , 0.0001          
-          1               , 0.00500468      , 0.0001          
-          2               , 0.00247269      , 0.0001          
-          3               , 0.000899035     , 0.0001          
+          0               , 0.00765473      , 0.0001          
+          1               , 0.00497838      , 0.0001          
+          2               , 0.0024697       , 0.0001          
+          3               , 0.00090216      , 0.0001          

TestCases/disc_adj_fsi/Airfoil_2d/configFlow.cfg

@@ -51,10 +51,9 @@ TIME_DISCRE_FLOW= EULER_IMPLICIT
 % Linear solvers ------------------------------------------------------- %
 LINEAR_SOLVER= FGMRES
 LINEAR_SOLVER_PREC= ILU
-LINEAR_SOLVER_ERROR= 1E-3
-LINEAR_SOLVER_ITER= 20
-LINEAR_SOLVER_SMOOTHER_RELAXATION= 0.7
-DISCADJ_LIN_SOLVER= SMOOTHER
+LINEAR_SOLVER_ERROR= 1E-4
+LINEAR_SOLVER_ITER= 50
+DISCADJ_LIN_SOLVER= FGMRES
 DISCADJ_LIN_PREC= ILU
 NEWTON_KRYLOV= YES
 QUASI_NEWTON_NUM_SAMPLES= 999
@@ -91,4 +90,4 @@ VOLUME_ADJ_FILENAME= adjoint_fluid
 SURFACE_FILENAME= surface_fluid
 SURFACE_ADJ_FILENAME= adjoint_surface_fluid
 SCREEN_WRT_FREQ_INNER= 10
-CONV_FILENAME= history
+CONV_FILENAME= history

TestCases/disc_adj_fsi/dyn_fsi/config.cfg

@@ -11,16 +11,15 @@ MESH_FILENAME= mesh.su2
 
 TIME_DOMAIN = YES
 TIME_ITER = 3
-TIME_STEP = 0.01
+TIME_STEP = 0.002
 UNST_ADJOINT_ITER= 10
 ITER_AVERAGE_OBJ= 3
 
 OBJECTIVE_FUNCTION= REFERENCE_NODE
 
-SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0]
+SCREEN_OUTPUT= TIME_ITER, OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], SENSITIVITY[1]
 
 WRT_ZONE_HIST=YES
 %WRT_ZONE_CONV=YES
 
-% Debug: Force Zero Grid Velocity = NO
 WRT_PERFORMANCE= NO

TestCases/disc_adj_fsi/dyn_fsi/configFEA.cfg

@@ -17,7 +17,7 @@ DESIGN_VARIABLE_FEA= YOUNG_MODULUS
 
 % Solid properties ----------------------------------------------------- %
 MATERIAL_MODEL= NEO_HOOKEAN
-ELASTICITY_MODULUS= 7E8
+ELASTICITY_MODULUS= 2E9
 POISSON_RATIO= 0.35
 MATERIAL_DENSITY= 2700.0
 
@@ -48,7 +48,6 @@ INCREMENTAL_LOAD= NO
 MESH_FILENAME= mesh.su2
 MESH_FORMAT= SU2
 
-RESTART_SOL= NO
 SOLUTION_FILENAME= solution_solid.dat
 SOLUTION_ADJ_FILENAME= adjoint_solid.dat
 OUTPUT_WRT_FREQ= 1

TestCases/disc_adj_fsi/dyn_fsi/configFlow.cfg

@@ -46,12 +46,15 @@ NUM_METHOD_GRAD= GREEN_GAUSS
 CFL_NUMBER= 1000.0
 
 % Flow numerics -------------------------------------------------------- %
-CONV_NUM_METHOD_FLOW= JST
+CONV_NUM_METHOD_FLOW= ROE
+MUSCL_FLOW= YES
+SLOPE_LIMITER_FLOW= NONE
 JST_SENSOR_COEFF= ( 0.5, 0.02 )
+ENTROPY_FIX_COEFF= 0.01
 TIME_DISCRE_FLOW= EULER_IMPLICIT
 
 % Linear solvers ------------------------------------------------------- %
-LINEAR_SOLVER= BCGSTAB
+LINEAR_SOLVER= FGMRES
 LINEAR_SOLVER_PREC= ILU
 LINEAR_SOLVER_ERROR= 1E-6
 LINEAR_SOLVER_ITER= 25
@@ -82,7 +85,6 @@ CONV_RESIDUAL_MINVAL= -7
 MESH_FILENAME= mesh.su2
 MESH_FORMAT= SU2
 %
-RESTART_SOL= NO
 SOLUTION_FILENAME= solution_fluid.dat
 SOLUTION_ADJ_FILENAME= adjoint_fluid.dat
 OUTPUT_WRT_FREQ= 1
@@ -95,5 +97,6 @@ SURFACE_FILENAME= surface_fluid
 SURFACE_ADJ_FILENAME= adjoint_surface_fluid
 SCREEN_WRT_FREQ_INNER= 10
 CONV_FILENAME= history
+SCREEN_OUTPUT= ITER, RMS_RES, LINSOL
 HISTORY_OUTPUT= ITER, RMS_RES, AERO_COEFF, TAVG_AERO_COEFF
 OUTPUT_FILES= NONE

TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref

@@ -1,9 +1,9 @@
 INDEX	GRAD
-0	-3.533163525398154e-03
-1	-1.879941699367452e-03
-2	-8.078436962531523e-04
-3	-2.794845540524826e-04
-4	-2.790760674146328e-04
-5	-8.045490756184532e-04
-6	-1.867166054373919e-03
-7	-3.498473613679959e-03
+0	-4.570869215186209e-04
+1	-2.401466750293265e-04
+2	-9.134698389661222e-05
+3	-1.628087055796442e-05
+4	-1.741052457131011e-05
+5	-9.462489787924037e-05
+6	-2.452466259888652e-04
+7	-4.635483632874935e-04

TestCases/flamelet/03_laminar_premixed_ch4_flame_cht_cfd/lam_prem_ch4_cht_cfd_master.cfg

@@ -43,13 +43,15 @@ SCREEN_WRT_FREQ_OUTER= 1
 
 WRT_VOLUME_OVERWRITE= YES
 
-WRT_ZONE_CONV= YES
+%WRT_ZONE_CONV= YES
 CONV_RESIDUAL_MINVAL= -20
 
 VOLUME_FILENAME= fluid
 
 CONV_FILENAME= history
 
+SCREEN_OUTPUT= OUTER_ITER, RMS_RES[0], RMS_RES[1]
+
 % --------------------------- Optimization Parameters --------------------------%
 
 MESH_FILENAME= mesh_unstructured_cht.su2

TestCases/flamelet/04_laminar_premixed_ch4_flame_cht_ad/lam_prem_ch4_cht_ad_master.cfg

@@ -44,10 +44,11 @@ OUTER_ITER = 1100
 OUTPUT_WRT_FREQ= 50
 SCREEN_WRT_FREQ_INNER= 1
 SCREEN_WRT_FREQ_OUTER= 1
+SCREEN_OUTPUT= OUTER_ITER, BGS_RES[0], BGS_RES[1]
 
 WRT_VOLUME_OVERWRITE= YES
 
-WRT_ZONE_CONV= YES
+%WRT_ZONE_CONV= YES
 CONV_RESIDUAL_MINVAL= -20
 
 VOLUME_FILENAME= fluid