diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 62ddab7e956d..3c8d2852fd38 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -600,8 +600,9 @@ class CConfig { Kappa_4th_Flow, /*!< \brief JST 4th order dissipation coefficient for flow equations. */ Kappa_2nd_Heat, /*!< \brief 2nd order dissipation coefficient for heat equation. */ Kappa_4th_Heat, /*!< \brief 4th order dissipation coefficient for heat equation. */ - Cent_Jac_Fix_Factor; /*!< \brief Multiply the dissipation contribution to the Jacobian of central schemes + Cent_Jac_Fix_Factor, /*!< \brief Multiply the dissipation contribution to the Jacobian of central schemes by this factor to make the global matrix more diagonal dominant. */ + Cent_Inc_Jac_Fix_Factor; /*!< \brief Multiply the dissipation contribution to the Jacobian of incompressible central schemes */ su2double Geo_Waterline_Location; /*!< \brief Location of the waterline. */ su2double Min_Beta_RoeTurkel, /*!< \brief Minimum value of Beta for the Roe-Turkel low Mach preconditioner. */ @@ -4625,6 +4626,12 @@ class CConfig { */ su2double GetCent_Jac_Fix_Factor(void) const { return Cent_Jac_Fix_Factor; } + /*! + * \brief Factor by which to multiply the dissipation contribution to Jacobians of incompressible central schemes. + * \return The factor. + */ + su2double GetCent_Inc_Jac_Fix_Factor(void) const { return Cent_Inc_Jac_Fix_Factor; } + /*! * \brief Get the kind of integration scheme (explicit or implicit) * for the adjoint flow equations. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 9c5cfca1e377..3b2dce6a0630 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1819,6 +1819,8 @@ void CConfig::SetConfig_Options() { addBoolOption("USE_ACCURATE_FLUX_JACOBIANS", Use_Accurate_Jacobians, false); /*!\brief CENTRAL_JACOBIAN_FIX_FACTOR \n DESCRIPTION: Improve the numerical properties (diagonal dominance) of the global Jacobian matrix, 3 to 4 is "optimum" (central schemes) \ingroup Config*/ addDoubleOption("CENTRAL_JACOBIAN_FIX_FACTOR", Cent_Jac_Fix_Factor, 4.0); + /*!\brief CENTRAL_JACOBIAN_FIX_FACTOR \n DESCRIPTION: Control numerical properties of the global Jacobian matrix using a multiplication factor for incompressible central schemes \ingroup Config*/ + addDoubleOption("CENTRAL_INC_JACOBIAN_FIX_FACTOR", Cent_Inc_Jac_Fix_Factor, 1.0); /*!\brief CONV_NUM_METHOD_ADJFLOW * \n DESCRIPTION: Convective numerical method for the adjoint solver. @@ -4646,12 +4648,6 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ } } - /*--- Rotating frame is not yet supported with the incompressible solver. ---*/ - - if ((Kind_Solver == INC_EULER || Kind_Solver == INC_NAVIER_STOKES || Kind_Solver == INC_RANS) && (Kind_GridMovement == ROTATING_FRAME)) { - SU2_MPI::Error("Support for rotating frame simulation not yet implemented for incompressible flows.", CURRENT_FUNCTION); - } - /*--- Assert that there are two markers being analyzed if the pressure drop objective function is selected. ---*/ diff --git a/SU2_CFD/include/numerics/flow/convection/centered.hpp b/SU2_CFD/include/numerics/flow/convection/centered.hpp index c7ee7551890b..d2e3116dad4b 100644 --- a/SU2_CFD/include/numerics/flow/convection/centered.hpp +++ b/SU2_CFD/include/numerics/flow/convection/centered.hpp @@ -252,6 +252,8 @@ class CCentLaxInc_Flow final : public CNumerics { variable_density, /*!< \brief Variable density incompressible flows. */ energy; /*!< \brief computation with the energy equation. */ + su2double fix_factor; /*!< \brief Fix factor for Jacobians. */ + su2double** Jacobian_i = nullptr; /*!< \brief The Jacobian w.r.t. point i after computation. */ su2double** Jacobian_j = nullptr; /*!< \brief The Jacobian w.r.t. point j after computation. */ @@ -312,6 +314,8 @@ class CCentJSTInc_Flow final : public CNumerics { variable_density, /*!< \brief Variable density incompressible flows. */ energy; /*!< \brief computation with the energy equation. */ + su2double fix_factor; /*!< \brief Fix factor for Jacobians. */ + su2double** Jacobian_i = nullptr; /*!< \brief The Jacobian w.r.t. point i after computation. */ su2double** Jacobian_j = nullptr; /*!< \brief The Jacobian w.r.t. point j after computation. */ diff --git a/SU2_CFD/src/numerics/flow/convection/centered.cpp b/SU2_CFD/src/numerics/flow/convection/centered.cpp index 5d9b36220df9..040c676d25ee 100644 --- a/SU2_CFD/src/numerics/flow/convection/centered.cpp +++ b/SU2_CFD/src/numerics/flow/convection/centered.cpp @@ -344,6 +344,7 @@ CCentLaxInc_Flow::CCentLaxInc_Flow(unsigned short val_nDim, unsigned short val_n variable_density = (config->GetKind_DensityModel() == VARIABLE); /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); + fix_factor = config->GetCent_Inc_Jac_Fix_Factor(); energy = config->GetEnergy_Equation(); /*--- Artificial dissipation part ---*/ @@ -531,8 +532,8 @@ CNumerics::ResidualType<> CCentLaxInc_Flow::ComputeResidual(const CConfig* confi for (jVar = 0; jVar < nVar; jVar++) { ProjFlux[iVar] += Precon[iVar][jVar]*Epsilon_0*Diff_V[jVar]*StretchingFactor*MeanLambda; if (implicit) { - Jacobian_i[iVar][jVar] += Precon[iVar][jVar]*Epsilon_0*StretchingFactor*MeanLambda; - Jacobian_j[iVar][jVar] -= Precon[iVar][jVar]*Epsilon_0*StretchingFactor*MeanLambda; + Jacobian_i[iVar][jVar] += fix_factor*Precon[iVar][jVar]*Epsilon_0*StretchingFactor*MeanLambda; + Jacobian_j[iVar][jVar] -= fix_factor*Precon[iVar][jVar]*Epsilon_0*StretchingFactor*MeanLambda; } } } @@ -564,6 +565,7 @@ CCentJSTInc_Flow::CCentJSTInc_Flow(unsigned short val_nDim, unsigned short val_n energy = config->GetEnergy_Equation(); /* A grid is defined as dynamic if there's rigid grid movement or grid deformation AND the problem is time domain */ dynamic_grid = config->GetDynamic_Grid(); + fix_factor = config->GetCent_Inc_Jac_Fix_Factor(); /*--- Artifical dissipation part ---*/ @@ -759,8 +761,8 @@ CNumerics::ResidualType<> CCentJSTInc_Flow::ComputeResidual(const CConfig* confi for (jVar = 0; jVar < nVar; jVar++) { ProjFlux[iVar] += Precon[iVar][jVar]*(Epsilon_2*Diff_V[jVar] - Epsilon_4*Diff_Lapl[jVar])*StretchingFactor*MeanLambda; if (implicit) { - Jacobian_i[iVar][jVar] += Precon[iVar][jVar]*(Epsilon_2 + Epsilon_4*su2double(Neighbor_i+1))*StretchingFactor*MeanLambda; - Jacobian_j[iVar][jVar] -= Precon[iVar][jVar]*(Epsilon_2 + Epsilon_4*su2double(Neighbor_j+1))*StretchingFactor*MeanLambda; + Jacobian_i[iVar][jVar] += fix_factor*Precon[iVar][jVar]*(Epsilon_2 + Epsilon_4*su2double(Neighbor_i+1))*StretchingFactor*MeanLambda; + Jacobian_j[iVar][jVar] -= fix_factor*Precon[iVar][jVar]*(Epsilon_2 + Epsilon_4*su2double(Neighbor_j+1))*StretchingFactor*MeanLambda; } } } diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index e221024c8ba8..48ca37687162 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -1527,9 +1527,9 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - /*--- Load the conservative variables ---*/ + /*--- Load the primitive variables ---*/ - numerics->SetConservative(nodes->GetSolution(iPoint), nullptr); + numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nullptr); /*--- Set incompressible density ---*/ @@ -3051,33 +3051,19 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver LinSysRes.AddBlock(iPoint, Residual); if (implicit) { - - SetPreconditioner(config, iPoint); - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_i[iVar][jVar] = Preconditioner[iVar][jVar]; - } - } - - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) { - if (config->GetTime_Marching() == DT_STEPPING_1ST) - Jacobian_i[iVar][jVar] *= Volume_nP1 / TimeStep; - if (config->GetTime_Marching() == DT_STEPPING_2ND) - Jacobian_i[iVar][jVar] *= (Volume_nP1*3.0)/(2.0*TimeStep); - } - } - - if (!energy) { - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar][nDim+1] = 0.0; - Jacobian_i[nDim+1][iVar] = 0.0; - } + for (iVar = 1; iVar < nVar; iVar++) { + if (config->GetTime_Marching() == DT_STEPPING_1ST) + Jacobian_i[iVar][iVar] = Volume_nP1 / TimeStep; + if (config->GetTime_Marching() == DT_STEPPING_2ND) + Jacobian_i[iVar][iVar] = (Volume_nP1*3.0)/(2.0*TimeStep); } + for (iDim = 0; iDim < nDim; iDim++) + Jacobian_i[iDim+1][iDim+1] = Density*Jacobian_i[iDim+1][iDim+1]; + if (energy) Jacobian_i[nDim+1][nDim+1] = Density*Cp*Jacobian_i[nDim+1][nDim+1]; Jacobian.AddBlock2Diag(iPoint, Jacobian_i); - } + } } @@ -3276,31 +3262,21 @@ void CIncEulerSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver to the dual time source term. ---*/ if (!energy) Residual[nDim+1] = 0.0; LinSysRes.AddBlock(iPoint, Residual); - if (implicit) { - SetPreconditioner(config, iPoint); - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) { - Jacobian_i[iVar][jVar] = Preconditioner[iVar][jVar]; - } - } - for (iVar = 0; iVar < nVar; iVar++) { - for (jVar = 0; jVar < nVar; jVar++) { - if (config->GetTime_Marching() == DT_STEPPING_1ST) - Jacobian_i[iVar][jVar] *= Volume_nP1 / TimeStep; - if (config->GetTime_Marching() == DT_STEPPING_2ND) - Jacobian_i[iVar][jVar] *= (Volume_nP1*3.0)/(2.0*TimeStep); - } + if (implicit) { + for (iVar = 1; iVar < nVar; iVar++) { + if (config->GetTime_Marching() == DT_STEPPING_1ST) + Jacobian_i[iVar][iVar] = Volume_nP1 / TimeStep; + if (config->GetTime_Marching() == DT_STEPPING_2ND) + Jacobian_i[iVar][iVar] = (Volume_nP1*3.0)/(2.0*TimeStep); } + for (iDim = 0; iDim < nDim; iDim++) + Jacobian_i[iDim+1][iDim+1] = Density*Jacobian_i[iDim+1][iDim+1]; + if (energy) Jacobian_i[nDim+1][nDim+1] = Density*Cp*Jacobian_i[nDim+1][nDim+1]; - if (!energy) { - for (iVar = 0; iVar < nVar; iVar++) { - Jacobian_i[iVar][nDim+1] = 0.0; - Jacobian_i[nDim+1][iVar] = 0.0; - } - } Jacobian.AddBlock2Diag(iPoint, Jacobian_i); } + } } diff --git a/SU2_PY/SU2/eval/gradients.py b/SU2_PY/SU2/eval/gradients.py index cabeb4cc112f..6141b20e23cd 100644 --- a/SU2_PY/SU2/eval/gradients.py +++ b/SU2_PY/SU2/eval/gradients.py @@ -838,6 +838,14 @@ def findiff( config, state=None ): name = su2io.expand_time(name,config) link.extend(name) + # files: restart solution for dual-time stepping first and second order + if 'RESTART_FILE_1' in files: + name = files['RESTART_FILE_1'] + pull.append(name) + if 'RESTART_FILE_2' in files: + name = files['RESTART_FILE_2'] + pull.append(name) + # files: target equivarea distribution if 'EQUIV_AREA' in special_cases and 'TARGET_EA' in files: pull.append(files['TARGET_EA']) diff --git a/SU2_PY/finite_differences.py b/SU2_PY/finite_differences.py index 77ed4ce302af..9f84af1475a6 100755 --- a/SU2_PY/finite_differences.py +++ b/SU2_PY/finite_differences.py @@ -77,6 +77,21 @@ def finite_differences( filename , # State state = SU2.io.State() state.find_files(config) + + # add restart files to state.FILES + if config.get('TIME_DOMAIN', 'NO') == 'YES' and config.get('RESTART_SOL', 'NO') == 'YES': + restart_name = config['RESTART_FILENAME'].split('.')[0] + restart_filename = restart_name + '_' + str(int(config['RESTART_ITER'])-1).zfill(5) + '.dat' + if not os.path.isfile(restart_filename): # throw, if restart files does not exist + sys.exit("Error: Restart file <" + restart_filename + "> not found.") + state['FILES']['RESTART_FILE_1'] = restart_filename + + # use only, if time integration is second order + if config.get('TIME_MARCHING', 'NO') == 'DUAL_TIME_STEPPING-2ND_ORDER': + restart_filename = restart_name + '_' + str(int(config['RESTART_ITER'])-2).zfill(5) + '.dat' + if not os.path.isfile(restart_filename): # throw, if restart files does not exist + sys.exit("Error: Restart file <" + restart_filename + "> not found.") + state['FILES']['RESTART_FILE_2'] =restart_filename # Finite Difference Gradients SU2.eval.gradients.findiff(config,state) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index de5032f126ef..f950ea8a6201 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -866,7 +866,7 @@ def main(): unst_inc_turb_naca0015_sa.cfg_dir = "unsteady/pitching_naca0015_rans_inc" unst_inc_turb_naca0015_sa.cfg_file = "config_incomp_turb_sa.cfg" unst_inc_turb_naca0015_sa.test_iter = 1 - unst_inc_turb_naca0015_sa.test_vals = [-2.991060, -6.866340, 1.476415, 0.419712] + unst_inc_turb_naca0015_sa.test_vals = [-3.004011, -6.876230, 1.487888, 0.421869] unst_inc_turb_naca0015_sa.su2_exec = "parallel_computation.py -f" unst_inc_turb_naca0015_sa.timeout = 1600 unst_inc_turb_naca0015_sa.tol = 0.00001 @@ -1195,8 +1195,8 @@ def main(): cht_incompressible_unsteady.cfg_dir = "coupled_cht/incomp_2d_unsteady" cht_incompressible_unsteady.cfg_file = "cht_2d_3cylinders.cfg" cht_incompressible_unsteady.test_iter = 2 - cht_incompressible_unsteady.test_vals = [-1.348104, -0.080392, -0.080391, -0.080391] #last 4 columns - cht_incompressible_unsteady.su2_exec = "SU2_CFD" + cht_incompressible_unsteady.test_vals = [-1.305471, -0.080372, -0.080376, -0.080372] #last 4 columns + cht_incompressible_unsteady.su2_exec = "mpirun -n 2 SU2_CFD" cht_incompressible_unsteady.timeout = 1600 cht_incompressible_unsteady.multizone = True cht_incompressible_unsteady.unsteady = True diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index a44c57b8dfe9..539b1c0250ab 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1001,7 +1001,7 @@ def main(): unst_inc_turb_naca0015_sa.cfg_dir = "unsteady/pitching_naca0015_rans_inc" unst_inc_turb_naca0015_sa.cfg_file = "config_incomp_turb_sa.cfg" unst_inc_turb_naca0015_sa.test_iter = 1 - unst_inc_turb_naca0015_sa.test_vals = [-2.994996, -6.869781, 1.434864, 0.416626] #last 4 columns + unst_inc_turb_naca0015_sa.test_vals = [-3.007635, -6.879789, 1.445300, 0.419281] #last 4 columns unst_inc_turb_naca0015_sa.su2_exec = "SU2_CFD" unst_inc_turb_naca0015_sa.timeout = 1600 unst_inc_turb_naca0015_sa.tol = 0.00001 @@ -1407,7 +1407,7 @@ def main(): cht_incompressible_unsteady.cfg_dir = "coupled_cht/incomp_2d_unsteady" cht_incompressible_unsteady.cfg_file = "cht_2d_3cylinders.cfg" cht_incompressible_unsteady.test_iter = 2 - cht_incompressible_unsteady.test_vals = [-1.348104, -0.080392, -0.080391, -0.080391] #last 4 columns + cht_incompressible_unsteady.test_vals = [-1.303588, -0.080377, -0.080380, -0.080377] #last 4 columns cht_incompressible_unsteady.su2_exec = "SU2_CFD" cht_incompressible_unsteady.timeout = 1600 cht_incompressible_unsteady.multizone = True