From 503205949f52b3e40aab99e915a668967c848ece Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 26 Sep 2025 22:03:08 -0700 Subject: [PATCH 1/9] newton krylov relaxation to help the linear solver converge --- Common/include/CConfig.hpp | 21 +++-- Common/include/option_structure.inl | 15 ++-- Common/src/CConfig.cpp | 88 +++++++++---------- Common/src/linear_algebra/CSysSolve.cpp | 2 +- .../integration/CNewtonIntegration.hpp | 7 +- .../numerics_simd/flow/convection/roe.hpp | 7 +- .../src/integration/CNewtonIntegration.cpp | 58 +++++++++--- SU2_CFD/src/solvers/CEulerSolver.cpp | 3 +- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 3 +- SU2_CFD/src/solvers/CNEMOEulerSolver.cpp | 5 +- SU2_CFD/src/solvers/CSolver.cpp | 5 ++ TestCases/vandv/rans/30p30n/config.cfg | 15 ++-- 12 files changed, 147 insertions(+), 82 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index d33eb626a1e1..4bb6998dd25e 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -431,7 +431,8 @@ class CConfig { bool UseVectorization; /*!< \brief Whether to use vectorized numerics schemes. */ bool NewtonKrylov; /*!< \brief Use a coupled Newton method to solve the flow equations. */ array NK_IntParam{{20, 3, 2}}; /*!< \brief Integer parameters for NK method. */ - array NK_DblParam{{-2.0, 0.1, -3.0, 1e-4}}; /*!< \brief Floating-point parameters for NK method. */ + array NK_DblParam{{-2.0, 0.1, -3.0, 1e-4, 1.0}}; /*!< \brief Floating-point parameters for NK method. */ + su2double NK_Relaxation = 1.0; unsigned short nMGLevels; /*!< \brief Number of multigrid levels (coarse levels). */ unsigned short nCFL; /*!< \brief Number of CFL, one for each multigrid level. */ @@ -1332,9 +1333,9 @@ class CConfig { template void addEnumListOption(const string name, unsigned short& input_size, Tfield*& option_field, const map& enum_map); - void addDoubleArrayOption(const string& name, const int size, su2double* option_field); + void addDoubleArrayOption(const string& name, int size, bool allow_fewer, su2double* option_field); - void addUShortArrayOption(const string& name, const int size, unsigned short* option_field); + void addUShortArrayOption(const string& name, int size, bool allow_fewer, unsigned short* option_field); void addDoubleListOption(const string& name, unsigned short & size, su2double * & option_field); @@ -4369,12 +4370,22 @@ class CConfig { /*! * \brief Get Newton-Krylov integer parameters. */ - array GetNewtonKrylovIntParam(void) const { return NK_IntParam; } + array GetNewtonKrylovIntParam() const { return NK_IntParam; } /*! * \brief Get Newton-Krylov floating-point parameters. */ - array GetNewtonKrylovDblParam(void) const { return NK_DblParam; } + array GetNewtonKrylovDblParam() const { return NK_DblParam; } + + /*! + * \brief Get the Newton-Krylov relaxation. + */ + su2double GetNewtonKrylovRelaxation() const { return NK_Relaxation; } + + /*! + * \brief Set the Newton-Krylov relaxation. + */ + void SetNewtonKrylovRelaxation(const su2double& relaxation) { NK_Relaxation = relaxation; } /*! * \brief Returns the Roe kappa (multipler of the dissipation term). diff --git a/Common/include/option_structure.inl b/Common/include/option_structure.inl index 863e3aadd730..5dba5c845b22 100644 --- a/Common/include/option_structure.inl +++ b/Common/include/option_structure.inl @@ -233,18 +233,19 @@ class COptionEnumList final : public COptionBase { template class COptionArray final : public COptionBase { - string name; // Identifier for the option - const int size; // Number of elements - Type* field; // Reference to the field + string name; // Identifier for the option + const int size; // Number of elements + const bool allow_fewer; // Allow smaller size + Type* field; // Reference to the field public: - COptionArray(string option_field_name, const int list_size, Type* option_field) - : name(option_field_name), size(list_size), field(option_field) {} + COptionArray(string option_field_name, const int list_size, const bool allow_fewer, Type* option_field) + : name(std::move(option_field_name)), size(list_size), allow_fewer(allow_fewer), field(option_field) {} string SetValue(const vector& option_value) override { COptionBase::SetValue(option_value); // Check that the size is correct - if (option_value.size() != (unsigned long)this->size) { + if ((option_value.size() != size_t(size) && allow_fewer) || option_value.size() > size_t(size)) { string newstring; newstring.append(this->name); newstring.append(": wrong number of arguments: "); @@ -258,7 +259,7 @@ class COptionArray final : public COptionBase { newstring.append(" found"); return newstring; } - for (int i = 0; i < this->size; i++) { + for (size_t i = 0; i < option_value.size(); i++) { istringstream is(option_value[i]); if (!(is >> field[i])) { return badValue(" array", this->name); diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 834b0c953d90..8872561e63d8 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -367,17 +367,17 @@ void CConfig::addEnumListOption(const string name, unsigned short& input_size, T option_map.insert( pair(name, val) ); } -void CConfig::addDoubleArrayOption(const string& name, const int size, su2double* option_field) { +void CConfig::addDoubleArrayOption(const string& name, const int size, const bool allow_fewer, su2double* option_field) { assert(option_map.find(name) == option_map.end()); all_options.insert(pair(name, true)); - COptionBase* val = new COptionArray(name, size, option_field); + COptionBase* val = new COptionArray(name, size, allow_fewer, option_field); option_map.insert(pair(name, val)); } -void CConfig::addUShortArrayOption(const string& name, const int size, unsigned short* option_field) { +void CConfig::addUShortArrayOption(const string& name, const int size, const bool allow_fewer, unsigned short* option_field) { assert(option_map.find(name) == option_map.end()); all_options.insert(pair(name, true)); - COptionBase* val = new COptionArray(name, size, option_field); + COptionBase* val = new COptionArray(name, size, allow_fewer, option_field); option_map.insert(pair(name, val)); } @@ -1167,7 +1167,7 @@ void CConfig::SetConfig_Options() { addBoolOption("BODY_FORCE", Body_Force, false); body_force[0] = 0.0; body_force[1] = 0.0; body_force[2] = 0.0; /* DESCRIPTION: Vector of body force values (BodyForce_X, BodyForce_Y, BodyForce_Z) */ - addDoubleArrayOption("BODY_FORCE_VECTOR", 3, body_force); + addDoubleArrayOption("BODY_FORCE_VECTOR", 3, false, body_force); /* DESCRIPTION: Apply a body force as a source term for periodic boundary conditions \n Options: NONE, PRESSURE_DROP, MASSFLOW \n DEFAULT: NONE \ingroup Config */ addEnumOption("KIND_STREAMWISE_PERIODIC", Kind_Streamwise_Periodic, Streamwise_Periodic_Map, ENUM_STREAMWISE_PERIODIC::NONE); @@ -1314,11 +1314,11 @@ void CConfig::SetConfig_Options() { /*--- Options related to temperature polynomial coefficients for fluid models. ---*/ /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ - addDoubleArrayOption("CP_POLYCOEFFS", N_POLY_COEFFS, cp_polycoeffs.data()); + addDoubleArrayOption("CP_POLYCOEFFS", N_POLY_COEFFS, false, cp_polycoeffs.data()); /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ - addDoubleArrayOption("MU_POLYCOEFFS", N_POLY_COEFFS, mu_polycoeffs.data()); + addDoubleArrayOption("MU_POLYCOEFFS", N_POLY_COEFFS, false, mu_polycoeffs.data()); /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ - addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, kt_polycoeffs.data()); + addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, false, kt_polycoeffs.data()); /*!\brief REYNOLDS_NUMBER \n DESCRIPTION: Reynolds number (non-dimensional, based on the free-stream values). Needed for viscous solvers. For incompressible solvers the Reynolds length will always be 1.0 \n DEFAULT: 0.0 \ingroup Config */ addDoubleOption("REYNOLDS_NUMBER", Reynolds, 0.0); @@ -1378,7 +1378,7 @@ void CConfig::SetConfig_Options() { addDoubleOption("INC_DENSITY_INIT", Inc_Density_Init, 1.2886); /*!\brief INC_VELOCITY_INIT \n DESCRIPTION: Initial velocity for incompressible flows (1.0,0,0 m/s by default) \ingroup Config*/ vel_init[0] = 1.0; vel_init[1] = 0.0; vel_init[2] = 0.0; - addDoubleArrayOption("INC_VELOCITY_INIT", 3, vel_init); + addDoubleArrayOption("INC_VELOCITY_INIT", 3, false, vel_init); /*!\brief INC_TEMPERATURE_INIT \n DESCRIPTION: Initial temperature for incompressible flows with the energy equation (288.15 K by default) \ingroup Config*/ addDoubleOption("INC_TEMPERATURE_INIT", Inc_Temperature_Init, 288.15); /*!\brief INC_NONDIM \n DESCRIPTION: Non-dimensionalization scheme for incompressible flows. \ingroup Config*/ @@ -1403,10 +1403,10 @@ void CConfig::SetConfig_Options() { /*!\brief FLAME_INIT_METHOD \n DESCRIPTION: Ignition method for flamelet solver \n DEFAULT: no ignition; cold flow only. */ addEnumOption("FLAME_INIT_METHOD", flamelet_ParsedOptions.ignition_method, Flamelet_Init_Map, FLAMELET_INIT_TYPE::NONE); /*!\brief FLAME_INIT \n DESCRIPTION: flame front initialization using the flamelet model \ingroup Config*/ - addDoubleArrayOption("FLAME_INIT", flamelet_ParsedOptions.flame_init.size(),flamelet_ParsedOptions.flame_init.begin()); + addDoubleArrayOption("FLAME_INIT", flamelet_ParsedOptions.flame_init.size(), false, flamelet_ParsedOptions.flame_init.begin()); /*!\brief SPARK_INIT \n DESCRIPTION: spark initialization using the flamelet model \ingroup Config*/ - addDoubleArrayOption("SPARK_INIT", flamelet_ParsedOptions.spark_init.size(), flamelet_ParsedOptions.spark_init.begin()); + addDoubleArrayOption("SPARK_INIT", flamelet_ParsedOptions.spark_init.size(), false, flamelet_ParsedOptions.spark_init.begin()); /*!\brief SPARK_REACTION_RATES \n DESCRIPTION: Net source term values applied to species within spark area during spark ignition. \ingroup Config*/ addDoubleListOption("SPARK_REACTION_RATES", flamelet_ParsedOptions.nspark, flamelet_ParsedOptions.spark_reaction_rates); @@ -1426,7 +1426,7 @@ void CConfig::SetConfig_Options() { vel_inf[0] = 1.0; vel_inf[1] = 0.0; vel_inf[2] = 0.0; /*!\brief FREESTREAM_VELOCITY\n DESCRIPTION: Free-stream velocity (m/s) */ - addDoubleArrayOption("FREESTREAM_VELOCITY", 3, vel_inf); + addDoubleArrayOption("FREESTREAM_VELOCITY", 3, false, vel_inf); /* DESCRIPTION: Free-stream viscosity (1.853E-5 Ns/m^2 (air), 0.798E-3 Ns/m^2 (water)) */ addDoubleOption("FREESTREAM_VISCOSITY", Viscosity_FreeStream, -1.0); /* DESCRIPTION: */ @@ -1623,7 +1623,7 @@ void CConfig::SetConfig_Options() { /*!\brief GILES_EXTRA_RELAXFACTOR \n DESCRIPTION: the 1st coeff the value of the under relaxation factor to apply to the shroud and hub, * the 2nd coefficient is the the percentage of span-wise height influenced by this extra under relaxation factor.*/ extrarelfac[0] = 0.1; extrarelfac[1] = 0.1; - addDoubleArrayOption("GILES_EXTRA_RELAXFACTOR", 2, extrarelfac); + addDoubleArrayOption("GILES_EXTRA_RELAXFACTOR", 2, false, extrarelfac); /*!\brief AVERAGE_PROCESS_TYPE \n DESCRIPTION: types of mixing process for averaging quantities at the boundaries. \n OPTIONS: see \link MixingProcess_Map \endlink \n DEFAULT: AREA_AVERAGE \ingroup Config*/ addEnumOption("MIXINGPLANE_INTERFACE_KIND", Kind_MixingPlaneInterface, MixingPlaneInterface_Map, NEAREST_SPAN); @@ -1637,13 +1637,13 @@ void CConfig::SetConfig_Options() { /*!\brief MIXEDOUT_COEFF \n DESCRIPTION: the 1st coeff is an under relaxation factor for the Newton method, * the 2nd coefficient is the tolerance for the Newton method, 3rd coefficient is the maximum number of * iteration for the Newton Method.*/ - addDoubleArrayOption("MIXEDOUT_COEFF", 3, mixedout_coeff); + addDoubleArrayOption("MIXEDOUT_COEFF", 3, false, mixedout_coeff); /*!\brief RAMP_MOTION_FRAME\n DESCRIPTION: option to ramp up or down the frame of motion velocity value*/ addBoolOption("RAMP_MOTION_FRAME", RampMotionFrame, false); rampMotionFrameCoeff[0] = 100.0; rampMotionFrameCoeff[1] = 1.0; rampMotionFrameCoeff[2] = 1000.0; /*!\brief RAMP_MOTION_FRAME_COEFF \n DESCRIPTION: the 1st coeff is the staring outlet value, * the 2nd coeff is the number of iterations for the update, 3rd is the number of total iteration till reaching the final outlet pressure value */ - addDoubleArrayOption("RAMP_MOTION_FRAME_COEFF", 3, rampMotionFrameCoeff); + addDoubleArrayOption("RAMP_MOTION_FRAME_COEFF", 3, false, rampMotionFrameCoeff); /* DESCRIPTION: AVERAGE_MACH_LIMIT is a limit value for average procedure based on the mass flux. */ addDoubleOption("AVERAGE_MACH_LIMIT", AverageMachLimit, 0.03); /*!\brief RAMP_OUTLET\n DESCRIPTION: option to ramp up or down the Giles outlet value*/ @@ -1651,7 +1651,7 @@ void CConfig::SetConfig_Options() { rampOutletCoeff[0] = 100000.0; rampOutletCoeff[1] = 1.0; rampOutletCoeff[2] = 1000.0; /*!\brief RAMP_OUTLET_COEFF \n DESCRIPTION: the 1st coeff is the staring outlet value, * the 2nd coeff is the number of iterations for the update, 3rd is the number of total iteration till reaching the final outlet pressure value */ - addDoubleArrayOption("RAMP_OUTLET_COEFF", 3, rampOutletCoeff); + addDoubleArrayOption("RAMP_OUTLET_COEFF", 3, false, rampOutletCoeff); /*!\brief MARKER_MIXINGPLANE \n DESCRIPTION: Identify the boundaries in which the mixing plane is applied. \ingroup Config*/ addStringListOption("MARKER_MIXINGPLANE_INTERFACE", nMarker_MixingPlaneInterface, Marker_MixingPlaneInterface); /*!\brief TURBULENT_MIXINGPLANE \n DESCRIPTION: Activate mixing plane also for turbulent quantities \ingroup Config*/ @@ -1724,14 +1724,14 @@ void CConfig::SetConfig_Options() { addBoolOption("ACTDISK_SU2_DEF", ActDisk_SU2_DEF, false); /* DESCRIPTION: Definition of the distortion rack (radial number of proves / circumferential density (degree) */ distortion[0] = 5.0; distortion[1] = 15.0; - addDoubleArrayOption("DISTORTION_RACK", 2, distortion); + addDoubleArrayOption("DISTORTION_RACK", 2, false, distortion); /* DESCRIPTION: Values of the box to impose a subsonic nacellle (mach, Pressure, Temperature) */ eng_val[0]=0.0; eng_val[1]=0.0; eng_val[2]=0.0; eng_val[3]=0.0; eng_val[4]=0.0; - addDoubleArrayOption("SUBSONIC_ENGINE_VALUES", 5, eng_val); + addDoubleArrayOption("SUBSONIC_ENGINE_VALUES", 5, false, eng_val); /* DESCRIPTION: Coordinates of the box to impose a subsonic nacellle cylinder (Xmin, Ymin, Zmin, Xmax, Ymax, Zmax, Radius) */ eng_cyl[0] = 0.0; eng_cyl[1] = 0.0; eng_cyl[2] = 0.0; eng_cyl[3] = 1E15; eng_cyl[4] = 1E15; eng_cyl[5] = 1E15; eng_cyl[6] = 1E15; - addDoubleArrayOption("SUBSONIC_ENGINE_CYL", 7, eng_cyl); + addDoubleArrayOption("SUBSONIC_ENGINE_CYL", 7, false, eng_cyl); /* DESCRIPTION: Engine exhaust boundary marker(s) Format: (nacelle exhaust marker, total nozzle temp, total nozzle pressure, ... )*/ addExhaustOption("MARKER_ENGINE_EXHAUST", nMarker_EngineExhaust, Marker_EngineExhaust, Exhaust_Temperature_Target, Exhaust_Pressure_Target); @@ -1753,7 +1753,7 @@ void CConfig::SetConfig_Options() { addBoolOption("SINE_LOAD", Sine_Load, false); sineload_coeff[0] = 0.0; sineload_coeff[1] = 0.0; sineload_coeff[2] = 0.0; /*!\brief SINE_LOAD_COEFF \n DESCRIPTION: the 1st coeff is the amplitude, the 2nd is the frequency, 3rd is the phase in radians */ - addDoubleArrayOption("SINE_LOAD_COEFF", 3, sineload_coeff); + addDoubleArrayOption("SINE_LOAD_COEFF", 3, false, sineload_coeff); /*!\brief RAMP_AND_RELEASE\n DESCRIPTION: release the load after applying the ramp*/ addBoolOption("RAMP_AND_RELEASE_LOAD", RampAndRelease, false); @@ -1776,9 +1776,9 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Use a Newton-Krylov method. */ addBoolOption("NEWTON_KRYLOV", NewtonKrylov, false); /* DESCRIPTION: Integer parameters {startup iters, precond iters, initial tolerance relaxation}. */ - addUShortArrayOption("NEWTON_KRYLOV_IPARAM", NK_IntParam.size(), NK_IntParam.data()); + addUShortArrayOption("NEWTON_KRYLOV_IPARAM", NK_IntParam.size(), true, NK_IntParam.data()); /* DESCRIPTION: Double parameters {startup residual drop, precond tolerance, full tolerance residual drop, findiff step}. */ - addDoubleArrayOption("NEWTON_KRYLOV_DPARAM", NK_DblParam.size(), NK_DblParam.data()); + addDoubleArrayOption("NEWTON_KRYLOV_DPARAM", NK_DblParam.size(), true, NK_DblParam.data()); /* DESCRIPTION: Number of samples for quasi-Newton methods. */ addUnsignedShortOption("QUASI_NEWTON_NUM_SAMPLES", nQuasiNewtonSamples, 0); @@ -1986,7 +1986,7 @@ void CConfig::SetConfig_Options() { addEnumOption("SLOPE_LIMITER_FLOW", Kind_SlopeLimit_Flow, Limiter_Map, LIMITER::VENKATAKRISHNAN); jst_coeff[0] = 0.5; jst_coeff[1] = 0.02; /*!\brief JST_SENSOR_COEFF \n DESCRIPTION: 2nd and 4th order artificial dissipation coefficients for the JST method \ingroup Config*/ - addDoubleArrayOption("JST_SENSOR_COEFF", 2, jst_coeff); + addDoubleArrayOption("JST_SENSOR_COEFF", 2, false, jst_coeff); /*!\brief LAX_SENSOR_COEFF \n DESCRIPTION: 1st order artificial dissipation coefficients for the Lax-Friedrichs method. \ingroup Config*/ addDoubleOption("LAX_SENSOR_COEFF", Kappa_1st_Flow, 0.15); /*!\brief USE_ACCURATE_FLUX_JACOBIANS \n DESCRIPTION: Use numerically computed Jacobians for AUSM+up(2) and SLAU(2) \ingroup Config*/ @@ -2007,7 +2007,7 @@ void CConfig::SetConfig_Options() { addEnumOption("SLOPE_LIMITER_ADJFLOW", Kind_SlopeLimit_AdjFlow, Limiter_Map, LIMITER::VENKATAKRISHNAN); jst_adj_coeff[0] = 0.5; jst_adj_coeff[1] = 0.02; /*!\brief ADJ_JST_SENSOR_COEFF \n DESCRIPTION: 2nd and 4th order artificial dissipation coefficients for the adjoint JST method. \ingroup Config*/ - addDoubleArrayOption("ADJ_JST_SENSOR_COEFF", 2, jst_adj_coeff); + addDoubleArrayOption("ADJ_JST_SENSOR_COEFF", 2, false, jst_adj_coeff); /*!\brief LAX_SENSOR_COEFF \n DESCRIPTION: 1st order artificial dissipation coefficients for the adjoint Lax-Friedrichs method. \ingroup Config*/ addDoubleOption("ADJ_LAX_SENSOR_COEFF", Kappa_1st_AdjFlow, 0.15); @@ -2073,7 +2073,7 @@ void CConfig::SetConfig_Options() { geo_loc[0] = 0.0; geo_loc[1] = 1.0; /* DESCRIPTION: Definition of the airfoil section */ - addDoubleArrayOption("GEO_BOUNDS", 2, geo_loc); + addDoubleArrayOption("GEO_BOUNDS", 2, false, geo_loc); /* DESCRIPTION: Identify the body to slice */ addEnumOption("GEO_DESCRIPTION", Geo_Description, Geo_Description_Map, WING); /* DESCRIPTION: Z location of the waterline */ @@ -2085,7 +2085,7 @@ void CConfig::SetConfig_Options() { nacelle_location[0] = 0.0; nacelle_location[1] = 0.0; nacelle_location[2] = 0.0; nacelle_location[3] = 0.0; nacelle_location[4] = 0.0; /* DESCRIPTION: Definition of the nacelle location (higlite coordinates, tilt angle, toe angle) */ - addDoubleArrayOption("GEO_NACELLE_LOCATION", 5, nacelle_location); + addDoubleArrayOption("GEO_NACELLE_LOCATION", 5, false, nacelle_location); /* DESCRIPTION: Output sectional forces for specified markers. */ addBoolOption("GEO_PLOT_STATIONS", Plot_Section_Forces, false); /* DESCRIPTION: Mode of the GDC code (analysis, or gradient) */ @@ -2163,11 +2163,11 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: List of the length of the RECTANGLE or BOX grid in the x,y,z directions. (default: (1.0,1.0,1.0) ). */ mesh_box_length[0] = 1.0; mesh_box_length[1] = 1.0; mesh_box_length[2] = 1.0; - addDoubleArrayOption("MESH_BOX_LENGTH", 3, mesh_box_length); + addDoubleArrayOption("MESH_BOX_LENGTH", 3, false, mesh_box_length); /* DESCRIPTION: List of the offset from 0.0 of the RECTANGLE or BOX grid in the x,y,z directions. (default: (0.0,0.0,0.0) ). */ mesh_box_offset[0] = 0.0; mesh_box_offset[1] = 0.0; mesh_box_offset[2] = 0.0; - addDoubleArrayOption("MESH_BOX_OFFSET", 3, mesh_box_offset); + addDoubleArrayOption("MESH_BOX_OFFSET", 3, false, mesh_box_offset); /* DESCRIPTION: Polynomial degree of the FEM solution for the RECTANGLE or BOX grid. (default: 1). */ addUnsignedShortOption("MESH_BOX_POLY_SOL_FEM", Mesh_Box_PSolFEM, 1); @@ -2257,21 +2257,21 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Mach number (non-dimensional, based on the mesh velocity and freestream vals.) */ addDoubleOption("MACH_MOTION", Mach_Motion, 0.0); /* DESCRIPTION: Coordinates of the rigid motion origin */ - addDoubleArrayOption("MOTION_ORIGIN", 3, Motion_Origin); + addDoubleArrayOption("MOTION_ORIGIN", 3, false, Motion_Origin); /* DESCRIPTION: Translational velocity vector (m/s) in the x, y, & z directions (RIGID_MOTION only) */ - addDoubleArrayOption("TRANSLATION_RATE", 3, Translation_Rate); + addDoubleArrayOption("TRANSLATION_RATE", 3, false, Translation_Rate); /* DESCRIPTION: Angular velocity vector (rad/s) about x, y, & z axes (RIGID_MOTION only) */ - addDoubleArrayOption("ROTATION_RATE", 3, Rotation_Rate); + addDoubleArrayOption("ROTATION_RATE", 3, false, Rotation_Rate); /* DESCRIPTION: Pitching angular freq. (rad/s) about x, y, & z axes (RIGID_MOTION only) */ - addDoubleArrayOption("PITCHING_OMEGA", 3, Pitching_Omega); + addDoubleArrayOption("PITCHING_OMEGA", 3, false, Pitching_Omega); /* DESCRIPTION: Pitching amplitude (degrees) about x, y, & z axes (RIGID_MOTION only) */ - addDoubleArrayOption("PITCHING_AMPL", 3, Pitching_Ampl); + addDoubleArrayOption("PITCHING_AMPL", 3, false, Pitching_Ampl); /* DESCRIPTION: Pitching phase offset (degrees) about x, y, & z axes (RIGID_MOTION only) */ - addDoubleArrayOption("PITCHING_PHASE", 3, Pitching_Phase); + addDoubleArrayOption("PITCHING_PHASE", 3, false, Pitching_Phase); /* DESCRIPTION: Plunging angular freq. (rad/s) in x, y, & z directions (RIGID_MOTION only) */ - addDoubleArrayOption("PLUNGING_OMEGA", 3, Plunging_Omega); + addDoubleArrayOption("PLUNGING_OMEGA", 3, false, Plunging_Omega); /* DESCRIPTION: Plunging amplitude (m) in x, y, & z directions (RIGID_MOTION only) */ - addDoubleArrayOption("PLUNGING_AMPL", 3, Plunging_Ampl); + addDoubleArrayOption("PLUNGING_AMPL", 3, false, Plunging_Ampl); /* DESCRIPTION: Coordinates of the rigid motion origin */ addDoubleListOption("SURFACE_MOTION_ORIGIN", nMarkerMotion_Origin, MarkerMotion_Origin); /* DESCRIPTION: Translational velocity vector (m/s) in the x, y, & z directions (DEFORMING only) */ @@ -2353,7 +2353,7 @@ void CConfig::SetConfig_Options() { addBoolOption("EQUIV_AREA", EquivArea, false); ea_lim[0] = 0.0; ea_lim[1] = 1.0; ea_lim[2] = 1.0; /* DESCRIPTION: Integration limits of the equivalent area ( xmin, xmax, Dist_NearField ) */ - addDoubleArrayOption("EA_INT_LIMIT", 3, ea_lim); + addDoubleArrayOption("EA_INT_LIMIT", 3, false, ea_lim); /* DESCRIPTION: Equivalent area scaling factor */ addDoubleOption("EA_SCALE_FACTOR", EA_ScaleFactor, 1.0); @@ -2402,7 +2402,7 @@ void CConfig::SetConfig_Options() { grid_fix[0] = -1E15; grid_fix[1] = -1E15; grid_fix[2] = -1E15; grid_fix[3] = 1E15; grid_fix[4] = 1E15; grid_fix[5] = 1E15; /* DESCRIPTION: Coordinates of the box where the grid will be deformed (Xmin, Ymin, Zmin, Xmax, Ymax, Zmax) */ - addDoubleArrayOption("HOLD_GRID_FIXED_COORD", 6, grid_fix); + addDoubleArrayOption("HOLD_GRID_FIXED_COORD", 6, false, grid_fix); /*!\par CONFIG_CATEGORY: Deformable mesh \ingroup Config*/ /*--- option related to deformable meshes ---*/ @@ -2521,7 +2521,7 @@ void CConfig::SetConfig_Options() { addDoubleOption("TOTAL_DV_PENALTY", DV_Penalty, 0); /*!\brief STRESS_PENALTY_PARAM\n DESCRIPTION: Maximum allowed stress and KS exponent for structural optimization \ingroup Config*/ - addDoubleArrayOption("STRESS_PENALTY_PARAM", 2, StressPenaltyParam.data()); + addDoubleArrayOption("STRESS_PENALTY_PARAM", 2, false, StressPenaltyParam.data()); /*!\brief REGIME_TYPE \n DESCRIPTION: Geometric condition \n OPTIONS: see \link Struct_Map \endlink \ingroup Config*/ addEnumOption("GEOMETRIC_CONDITIONS", Kind_Struct_Solver, Struct_Map, STRUCT_DEFORMATION::SMALL); @@ -2567,7 +2567,7 @@ void CConfig::SetConfig_Options() { inc_crit[0] = 0.0; inc_crit[1] = 0.0; inc_crit[2] = 0.0; /* DESCRIPTION: Definition of the UTOL RTOL ETOL*/ - addDoubleArrayOption("INCREMENTAL_CRITERIA", 3, inc_crit); + addDoubleArrayOption("INCREMENTAL_CRITERIA", 3, false, inc_crit); /* DESCRIPTION: Use of predictor */ addBoolOption("PREDICTOR", Predictor, false); @@ -2706,10 +2706,10 @@ void CConfig::SetConfig_Options() { addDoubleOption("HEAT_SOURCE_ROTATION_Z", Heat_Source_Rot_Z, 0.0); /* DESCRIPTION: Position of heat source center (Heat_Source_Center_X, Heat_Source_Center_Y, Heat_Source_Center_Z) */ hs_center[0] = 0.0; hs_center[1] = 0.0; hs_center[2] = 0.0; - addDoubleArrayOption("HEAT_SOURCE_CENTER", 3, hs_center); + addDoubleArrayOption("HEAT_SOURCE_CENTER", 3, false, hs_center); /* DESCRIPTION: Vector of heat source radii (Heat_Source_Axes_A, Heat_Source_Axes_B, Heat_Source_Axes_C) */ hs_axes[0] = 1.0; hs_axes[1] = 1.0; hs_axes[2] = 1.0; - addDoubleArrayOption("HEAT_SOURCE_AXES", 3, hs_axes); + addDoubleArrayOption("HEAT_SOURCE_AXES", 3, false, hs_axes); /*!\brief MARKER_EMISSIVITY DESCRIPTION: Wall emissivity of the marker for radiation purposes \n * Format: ( marker, emissivity of the marker, ... ) \ingroup Config */ @@ -2769,7 +2769,7 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Axis information for the spherical and cylindrical coord system */ ffd_axis[0] = 0.0; ffd_axis[1] = 0.0; ffd_axis[2] =0.0; - addDoubleArrayOption("FFD_AXIS", 3, ffd_axis); + addDoubleArrayOption("FFD_AXIS", 3, false, ffd_axis); /* DESCRIPTION: Number of total iterations in the FFD point inversion */ addUnsignedShortOption("FFD_ITERATIONS", nFFD_Iter, 500); @@ -2809,7 +2809,7 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Order of the BSplines for BSpline Blending function */ ffd_coeff[0] = 2; ffd_coeff[1] = 2; ffd_coeff[2] = 2; - addDoubleArrayOption("FFD_BSPLINE_ORDER", 3, ffd_coeff); + addDoubleArrayOption("FFD_BSPLINE_ORDER", 3, false, ffd_coeff); /*--- Options for the automatic differentiation methods ---*/ /*!\par CONFIG_CATEGORY: Automatic Differentation options\ingroup Config*/ diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp index fb0049e21d58..1a305c123dd3 100644 --- a/Common/src/linear_algebra/CSysSolve.cpp +++ b/Common/src/linear_algebra/CSysSolve.cpp @@ -47,7 +47,7 @@ constexpr T linSolEpsilon() { } template <> constexpr float linSolEpsilon() { - return 1e-12; + return 1e-14; } } // namespace diff --git a/SU2_CFD/include/integration/CNewtonIntegration.hpp b/SU2_CFD/include/integration/CNewtonIntegration.hpp index 450dc7934b26..f11d4dadff47 100644 --- a/SU2_CFD/include/integration/CNewtonIntegration.hpp +++ b/SU2_CFD/include/integration/CNewtonIntegration.hpp @@ -67,8 +67,10 @@ class CNewtonIntegration final : public CIntegration { enum class ResEvalType {EXPLICIT, DEFAULT}; bool setup = false; + bool autoRelaxation = false; Scalar finDiffStepND = 0.0; Scalar finDiffStep = 0.0; /*!< \brief Based on RMS(solution), used in matrix-free products. */ + Scalar nkRelaxation = 1.0; unsigned long omp_chunk_size; /*!< \brief Chunk size used in light point loops. */ /*--- Number of iterations and tolerance for the linear preconditioner, @@ -81,7 +83,7 @@ class CNewtonIntegration final : public CIntegration { * criteria are zero, or the solver does not provide a linear * preconditioner, there is no startup phase. ---*/ bool startupPeriod = false; - unsigned short startupIters = 0; + unsigned short startupIters = 0, iter = 0; su2double startupResidual = 0.0; su2double firstResidual = -20.0; @@ -96,8 +98,9 @@ class CNewtonIntegration final : public CIntegration { CNumerics*** numerics = nullptr; /*--- Residual and linear solver. ---*/ - CSysVector LinSysRes; + CSysVector LinSysRes, LinSysResRelax; CSysSolve LinSolver; + const CSysVector* LinSysRes0 = nullptr; /*--- If possible the solution vector of the solver is re-used, otherwise this temporary is used. ---*/ CSysVector LinSysSol; diff --git a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp index 7309c0f3f678..d46cacea3766 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp @@ -96,6 +96,7 @@ class CRoeBase : public Base { AD::StartPreacc(); const bool implicit = (config.GetKind_TimeIntScheme() == EULER_IMPLICIT); + const auto nkRelax = config.GetNewtonKrylovRelaxation(); const auto& solution = static_cast(solution_); const auto iPoint = geometry.edges->GetNode(iEdge,0); @@ -118,8 +119,12 @@ class CRoeBase : public Base { V1st.i.all = gatherVariables(iPoint, solution.GetPrimitive()); V1st.j.all = gatherVariables(jPoint, solution.GetPrimitive()); + VectorDbl mod_vector_ij; + for (size_t iDim = 0; iDim < nDim; ++iDim) { + mod_vector_ij(iDim) = nkRelax * vector_ij(iDim); + } auto V = reconstructPrimitives >( - iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, vector_ij, solution); + iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, mod_vector_ij, solution); /*--- Compute conservative variables. ---*/ diff --git a/SU2_CFD/src/integration/CNewtonIntegration.cpp b/SU2_CFD/src/integration/CNewtonIntegration.cpp index daaa7fcd64fc..455e239d3039 100644 --- a/SU2_CFD/src/integration/CNewtonIntegration.cpp +++ b/SU2_CFD/src/integration/CNewtonIntegration.cpp @@ -66,13 +66,18 @@ void CNewtonIntegration::Setup() { auto iparam = config->GetNewtonKrylovIntParam(); auto dparam = config->GetNewtonKrylovDblParam(); - startupIters = iparam[0]; + startupIters = iter = iparam[0]; startupResidual = dparam[0]; precondIters = iparam[1]; precondTol = dparam[1]; tolRelaxFactor = iparam[2]; fullTolResidual = dparam[2]; finDiffStepND = SU2_TYPE::GetValue(dparam[3]); + nkRelaxation = fmin(SU2_TYPE::GetValue(dparam[4]), 1); + if (nkRelaxation < 0) { + autoRelaxation = true; + nkRelaxation = 1; + } const auto nVar = solvers[FLOW_SOL]->GetnVar(); const auto nPoint = geometry->GetnPoint(); @@ -83,6 +88,9 @@ void CNewtonIntegration::Setup() { LinSolver.SetxIsZero(true); LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); + if (autoRelaxation || nkRelaxation < 1) { + LinSysResRelax.Initialize(nPoint, nPointDomain, nVar, 0.0); + } if (!std::is_same::value) { LinSysSol.Initialize(nPoint, nPointDomain, nVar, nullptr); @@ -175,6 +183,17 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver ** if (!setup) { Setup(); setup = true; } + // Ramp from 1st to 2nd order during the startup. + su2double baseNkRelaxation = 1; + if (startupPeriod) { + baseNkRelaxation = su2double(startupIters - iter) / startupIters; + } + config->SetNewtonKrylovRelaxation(baseNkRelaxation); + + // When using NK relaxation (not fully 2nd order Jacobian products) we need an additional + // residual evaluation that is used as the reference for finite differences. + LinSysRes0 = (!startupPeriod && nkRelaxation < 1) ? &LinSysResRelax : &LinSysRes; + SU2_OMP_PARALLEL_(if(solvers[FLOW_SOL]->GetHasHybridParallel())) { /*--- Save the current solution to be able to perturb it. ---*/ @@ -191,10 +210,20 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver ** if (preconditioner) preconditioner->Build(); - SU2_OMP_FOR_STAT(omp_chunk_size) - for (auto i = 0ul; i < LinSysRes.GetNElmDomain(); ++i) - LinSysRes[i] = SU2_TYPE::GetValue(solvers[FLOW_SOL]->LinSysRes[i]); - END_SU2_OMP_FOR + auto CopyLinSysRes = [&](int sign, auto& dst) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto i = 0ul; i < dst.GetNElmDomain(); ++i) + dst[i] = sign * SU2_TYPE::GetValue(solvers[FLOW_SOL]->LinSysRes[i]); + END_SU2_OMP_FOR + }; + CopyLinSysRes(1, LinSysRes); + + if (!startupPeriod && nkRelaxation < 1) { + SU2_OMP_SAFE_GLOBAL_ACCESS(config->SetNewtonKrylovRelaxation(nkRelaxation);) + ComputeResiduals(ResEvalType::EXPLICIT); + // Here the sign was not flipped by PrepareImplicitIteration. + CopyLinSysRes(-1, LinSysResRelax); + } su2double residual = 0.0; for (auto iVar = 0ul; iVar < LinSysRes.GetNVar(); ++iVar) @@ -207,10 +236,10 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver ** if (startupPeriod) { BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { firstResidual = max(firstResidual, residual); - if (startupIters) startupIters -= 1; + if (iter) iter -= 1; } END_SU2_OMP_SAFE_GLOBAL_ACCESS - endStartup = (startupIters == 0) && (residual - firstResidual < startupResidual); + endStartup = (iter == 0) && (residual - firstResidual < startupResidual); } /*--- The NK solves are expensive, the tolerance is relaxed while the residuals are high. ---*/ @@ -232,8 +261,7 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver ** if (startupPeriod) { iter = Preconditioner_impl(LinSysRes, linSysSol, iter, eps); - } - else { + } else { ComputeFinDiffStep(); eps *= toleranceFactor; @@ -247,6 +275,14 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver ** BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS { solvers[FLOW_SOL]->SetIterLinSolver(iter); solvers[FLOW_SOL]->SetResLinSolver(eps); + + if (!startupPeriod && autoRelaxation) { + if (eps > config->GetLinear_Solver_Error()) { + nkRelaxation *= 0.9; + } else if (iter < config->GetLinear_Solver_Iter()) { + nkRelaxation = fmin(nkRelaxation * 1.05, 1); + } + } } END_SU2_OMP_SAFE_GLOBAL_ACCESS @@ -303,11 +339,11 @@ void CNewtonIntegration::MatrixFreeProduct(const CSysVector& u, CSysVect su2double delta = (geometry->nodes->GetVolume(iPoint) + geometry->nodes->GetPeriodicVolume(iPoint)) / max(EPS, solvers[FLOW_SOL]->GetNodes()->GetDelta_Time(iPoint)); SU2_OMP_SIMD - for (auto iVar = 0ul; iVar < LinSysRes.GetNVar(); ++iVar) { + for (auto iVar = 0ul; iVar < LinSysRes0->GetNVar(); ++iVar) { Scalar perturbRes = SU2_TYPE::GetValue(solvers[FLOW_SOL]->LinSysRes(iPoint,iVar)); /*--- The global residual had its sign flipped, so we add to get the difference. ---*/ - v(iPoint,iVar) = (perturbRes + LinSysRes(iPoint,iVar)) * factor; + v(iPoint,iVar) = (perturbRes + (*LinSysRes0)(iPoint,iVar)) * factor; /*--- Pseudotime term of the true Jacobian. ---*/ v(iPoint,iVar) += SU2_TYPE::GetValue(delta) * u(iPoint,iVar); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 5368c287fc6b..8453e0fd7c84 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1796,6 +1796,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain } const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); + const su2double nkRelax = config->GetNewtonKrylovRelaxation(); const bool roe_turkel = (config->GetKind_Upwind_Flow() == UPWIND::TURKEL); const auto kind_dissipation = config->GetKind_RoeLowDiss(); @@ -1875,7 +1876,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain su2double Vector_ij[MAXNDIM] = {0.0}; for (iDim = 0; iDim < nDim; iDim++) { - Vector_ij[iDim] = 0.5*(Coord_j[iDim] - Coord_i[iDim]); + Vector_ij[iDim] = nkRelax * 0.5 * (Coord_j[iDim] - Coord_i[iDim]); } auto Gradient_i = nodes->GetGradient_Reconstruction(iPoint); diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 0910edbc5ebc..f4459a81bb89 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -1231,6 +1231,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE); const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); const bool bounded_scalar = config->GetBounded_Scalar(); + const su2double nkRelax = config->GetNewtonKrylovRelaxation(); /*--- For hybrid parallel AD, pause preaccumulation if there is shared reading of * variables, otherwise switch to the faster adjoint evaluation mode. ---*/ @@ -1271,7 +1272,7 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont su2double Vector_ij[MAXNDIM] = {0.0}; for (iDim = 0; iDim < nDim; iDim++) { - Vector_ij[iDim] = 0.5*(Coord_j[iDim] - Coord_i[iDim]); + Vector_ij[iDim] = nkRelax * 0.5 * (Coord_j[iDim] - Coord_i[iDim]); } auto Gradient_i = nodes->GetGradient_Reconstruction(iPoint); diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 27a77f5cb7e2..84c075d40d31 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -467,6 +467,7 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con const bool muscl = (config->GetMUSCL_Flow() && (iMesh == MESH_0)); const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE); const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); + const su2double nkRelax = config->GetNewtonKrylovRelaxation(); /*--- Non-physical counter. ---*/ unsigned long counter_local = 0; @@ -521,7 +522,7 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con /*--- High order reconstruction using MUSCL strategy ---*/ su2double Vector_ij[MAXNDIM] = {0.0}; for (auto iDim = 0ul; iDim < nDim; iDim++) { - Vector_ij[iDim] = 0.5*(Coord_j[iDim] - Coord_i[iDim]); + Vector_ij[iDim] = nkRelax * 0.5 * (Coord_j[iDim] - Coord_i[iDim]); } /*--- Retrieve gradient information ---*/ @@ -933,7 +934,7 @@ void CNEMOEulerSolver::ComputeUnderRelaxationFactor(const CConfig *config) { /* Loop over the solution update given by relaxing the linear system for this nonlinear iteration. */ - const su2double allowableRatio = 0.2; + const su2double allowableRatio = config->GetMaxUpdateFractionFlow(); SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 271859c7244f..e28e7d281978 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1779,19 +1779,24 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry, /* Sum the RMS residuals for all equations. */ New_Func = 0.0; + unsigned short totalVars = 0; for (unsigned short iVar = 0; iVar < solverFlow->GetnVar(); iVar++) { New_Func += log10(solverFlow->GetRes_RMS(iVar)); + ++totalVars; } if ((iMesh == MESH_0) && solverTurb) { for (unsigned short iVar = 0; iVar < solverTurb->GetnVar(); iVar++) { New_Func += log10(solverTurb->GetRes_RMS(iVar)); + ++totalVars; } } if ((iMesh == MESH_0) && solverSpecies) { for (unsigned short iVar = 0; iVar < solverSpecies->GetnVar(); iVar++) { New_Func += log10(solverSpecies->GetRes_RMS(iVar)); + ++totalVars; } } + New_Func /= totalVars; /* Compute the difference in the nonlinear residuals between the current and previous iterations, taking care with very low initial diff --git a/TestCases/vandv/rans/30p30n/config.cfg b/TestCases/vandv/rans/30p30n/config.cfg index da6749b88b36..385aede89fcb 100644 --- a/TestCases/vandv/rans/30p30n/config.cfg +++ b/TestCases/vandv/rans/30p30n/config.cfg @@ -60,23 +60,24 @@ MUSCL_TURB= NO TIME_DISCRE_FLOW= EULER_IMPLICIT TIME_DISCRE_TURB= EULER_IMPLICIT % -CFL_NUMBER= 40 +CFL_NUMBER= 10 CFL_REDUCTION_TURB= 1.0 -CFL_ADAPT= NO +CFL_ADAPT= YES +CFL_ADAPT_PARAM= ( 1, 1.05, 10, 150 ) % LINEAR_SOLVER= FGMRES LINEAR_SOLVER_PREC= ILU LINEAR_SOLVER_ERROR= 0.2 -LINEAR_SOLVER_ITER= 4 +LINEAR_SOLVER_ITER= 6 % MGLEVEL= 0 NEWTON_KRYLOV= YES -NEWTON_KRYLOV_IPARAM= ( 0, 0, 1 ) % n0, np, ft -NEWTON_KRYLOV_DPARAM= ( 0.0, 1e-20, -3, 1e-5 ) % r0, tp, rf, e +NEWTON_KRYLOV_IPARAM= ( 500, 0, 1 ) % n0, np, ft +NEWTON_KRYLOV_DPARAM= ( 0.0, 1e-20, -3, 1e-5, -1 ) % r0, tp, rf, e % % ------------------------ CONVERGENCE CRITERIA ------------------------- % % -ITER= 20 +ITER= 20000 CONV_RESIDUAL_MINVAL= -11.5 % % --------------------------- INPUT / OUTPUT ---------------------------- % @@ -89,5 +90,5 @@ SCREEN_WRT_FREQ_INNER= 1 HISTORY_OUTPUT= ( ITER, RMS_RES, AERO_COEFF ) SCREEN_OUTPUT= ( INNER_ITER, RMS_DENSITY, RMS_MOMENTUM-X, RMS_MOMENTUM-Y,\ RMS_ENERGY, RMS_NU_TILDE, DRAG, LIFT, MOMENT_Z,\ - LINSOL_RESIDUAL ) + LINSOL_RESIDUAL, LINSOL_ITER, AVG_CFL ) From e1befed52fa22b96d1c1290df05a32e820993810 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 27 Sep 2025 11:13:45 -0700 Subject: [PATCH 2/9] add badge --- README.md | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/README.md b/README.md index 9d17c535f1d4..7abff1642d9b 100644 --- a/README.md +++ b/README.md @@ -7,12 +7,10 @@ Computational analysis tools have revolutionized the way we design engineering systems, but most established codes are proprietary, unavailable, or prohibitively expensive for many users. The SU2 team is changing this, making multiphysics analysis and design optimization freely available as open-source software and involving everyone in its creation and development. -For an overview of the technical details in SU2, please see the following AIAA Journal article: - -"SU2: An open-source suite for multiphysics simulation and design," AIAA Journal, 54(3):828-846, 2016. - Please note that this project is released with a [Contributor Code of Conduct](CODE_OF_CONDUCT.md). By participating in this project you agree to abide by its terms. +If you use the latest version of SU2, please cite [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.17209057.svg)](https://doi.org/10.5281/zenodo.17209057) to acknowledge the active contributors. + Continuous Integration:
[![Regression Testing](https://github.com/su2code/SU2/workflows/Regression%20Testing/badge.svg?branch=develop)](https://github.com/su2code/SU2/actions) [![Release](https://github.com/su2code/SU2/workflows/Release%20Management/badge.svg?branch=develop)](https://github.com/su2code/SU2/actions) From 7059a78eeb7ac1349a6d28f7cd38ff1e1cad6590 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 27 Sep 2025 19:20:25 -0700 Subject: [PATCH 3/9] update docs --- Common/src/CConfig.cpp | 2 +- SU2_CFD/src/integration/CNewtonIntegration.cpp | 2 +- TestCases/euler/ramp/inv_ramp.cfg | 14 +++++++------- TestCases/rans/oneram6/turb_ONERAM6_nk.cfg | 17 ++++++++++++----- TestCases/vandv/rans/30p30n/config.cfg | 12 ++++++++---- config_template.cfg | 7 ++++--- 6 files changed, 33 insertions(+), 21 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 8872561e63d8..1d9700d6b689 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1777,7 +1777,7 @@ void CConfig::SetConfig_Options() { addBoolOption("NEWTON_KRYLOV", NewtonKrylov, false); /* DESCRIPTION: Integer parameters {startup iters, precond iters, initial tolerance relaxation}. */ addUShortArrayOption("NEWTON_KRYLOV_IPARAM", NK_IntParam.size(), true, NK_IntParam.data()); - /* DESCRIPTION: Double parameters {startup residual drop, precond tolerance, full tolerance residual drop, findiff step}. */ + /* DESCRIPTION: Double parameters {startup residual drop, precond tolerance, full tolerance residual drop, findiff step, NK relaxation}. */ addDoubleArrayOption("NEWTON_KRYLOV_DPARAM", NK_DblParam.size(), true, NK_DblParam.data()); /* DESCRIPTION: Number of samples for quasi-Newton methods. */ diff --git a/SU2_CFD/src/integration/CNewtonIntegration.cpp b/SU2_CFD/src/integration/CNewtonIntegration.cpp index 455e239d3039..a3286fd7c012 100644 --- a/SU2_CFD/src/integration/CNewtonIntegration.cpp +++ b/SU2_CFD/src/integration/CNewtonIntegration.cpp @@ -185,7 +185,7 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver ** // Ramp from 1st to 2nd order during the startup. su2double baseNkRelaxation = 1; - if (startupPeriod) { + if (startupPeriod && startupIters > 0 && !config->GetRestart()) { baseNkRelaxation = su2double(startupIters - iter) / startupIters; } config->SetNewtonKrylovRelaxation(baseNkRelaxation); diff --git a/TestCases/euler/ramp/inv_ramp.cfg b/TestCases/euler/ramp/inv_ramp.cfg index a5d3fb5b77cc..ee9aa387bf5a 100644 --- a/TestCases/euler/ramp/inv_ramp.cfg +++ b/TestCases/euler/ramp/inv_ramp.cfg @@ -45,21 +45,21 @@ NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES CFL_NUMBER= 1.0 CFL_ADAPT= YES CFL_ADAPT_PARAM= ( 0.9, 1.1, 1.0, 100.0, 0.8) -ITER= 500 -LINEAR_SOLVER= FGMRES +ITER= 500 +LINEAR_SOLVER= FGMRES LINEAR_SOLVER_ERROR= 1E-6 LINEAR_SOLVER_ITER= 3 NEWTON_KRYLOV= YES NEWTON_KRYLOV_IPARAM= (0, 1, 1) % n0, np, ft -NEWTON_KRYLOV_DPARAM= (0.0, 1e-20, -2.0, 1e-5) % r0, tp, rf, e +NEWTON_KRYLOV_DPARAM= (0.0, 1e-20, -2.0, 1e-5, 1.0) % r0, tp, rf, e, nkr MGLEVEL= 0 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% % CONV_NUM_METHOD_FLOW= AUSMPLUSUP MUSCL_FLOW= YES -SLOPE_LIMITER_FLOW= NISHIKAWA_R5 -VENKAT_LIMITER_COEFF= 0.5 +SLOPE_LIMITER_FLOW= NISHIKAWA_R5 +VENKAT_LIMITER_COEFF= 0.5 TIME_DISCRE_FLOW= EULER_IMPLICIT % --------------------------- CONVERGENCE PARAMETERS --------------------------% @@ -71,8 +71,8 @@ CONV_CAUCHY_EPS= 1E-10 % ------------------------- INPUT/OUTPUT INFORMATION --------------------------% % -MESH_FILENAME= ramp_unst.su2 -MESH_FORMAT= SU2 +MESH_FILENAME= ramp_unst.su2 +MESH_FORMAT= SU2 MESH_OUT_FILENAME= mesh_out SOLUTION_FILENAME= restart_flow SOLUTION_ADJ_FILENAME= solution_adj diff --git a/TestCases/rans/oneram6/turb_ONERAM6_nk.cfg b/TestCases/rans/oneram6/turb_ONERAM6_nk.cfg index 3e4c9caac60b..dd1f7fdfdbff 100644 --- a/TestCases/rans/oneram6/turb_ONERAM6_nk.cfg +++ b/TestCases/rans/oneram6/turb_ONERAM6_nk.cfg @@ -18,8 +18,8 @@ NEWTON_KRYLOV= YES % Iterations and tolerance for the Krylov part, it is important not to % "over solve", tolerance should be as high as possible. -LINEAR_SOLVER_ITER= 5 -LINEAR_SOLVER_ERROR= 0.25 +LINEAR_SOLVER_ITER= 8 +LINEAR_SOLVER_ERROR= 0.2 % For "n0" iterations or "r0" residual reduction, the normal quasi-Newton iterations % are used. Then, they become the preconditioner for the NK iterations with "np" linear @@ -27,13 +27,20 @@ LINEAR_SOLVER_ERROR= 0.25 % used directly (this may be enough for unsteady). % The tolerance for NK iterations is initially relaxed by factor "ft", and reaches % LINEAR_SOLVER_ERROR after "rf" residual reduction (additional to "r0"). +% If "ft"=1, "rf" does not matter. % The Jacobian-free products are based on finite differences with step "e". -NEWTON_KRYLOV_IPARAM= (10, 3, 2) % n0, np, ft -NEWTON_KRYLOV_DPARAM= (1.0, 0.1, -6.0, 1e-5) % r0, tp, rf, e +% "rnk" makes the linear system easier to solve by relaxing MUSCL (just for the NK +% matrix-free system). 0 is the same as not using NK, 1 is full MUSCL, a negative +% values means "automatic", the solver changes the value so that the linear system +% can be solved. +% If "n0">0, and RESTART_SOL=NO, the NK relaxation is applied to the residual (not +% just to the NK system) to ramp from MUSCL=OFF to ON during "n0" iterations. +NEWTON_KRYLOV_IPARAM= (200, 0, 1) % n0, np, ft +NEWTON_KRYLOV_DPARAM= (0, 0, 0, 1e-5, -1) % r0, tp, rf, e, rnk CFL_ADAPT= YES % it's needed CFL_NUMBER= 10 -CFL_ADAPT_PARAM= ( 0.8, 1.1, 5, 1000 ) % no point using NK with low CFL values +CFL_ADAPT_PARAM= ( 1, 1.05, 10, 400 ) % no point using NK with low CFL values % It is important (more than usual) to have similar magnitude variables REF_DIMENSIONALIZATION= FREESTREAM_VEL_EQ_MACH diff --git a/TestCases/vandv/rans/30p30n/config.cfg b/TestCases/vandv/rans/30p30n/config.cfg index 385aede89fcb..f4b9fa384a37 100644 --- a/TestCases/vandv/rans/30p30n/config.cfg +++ b/TestCases/vandv/rans/30p30n/config.cfg @@ -63,22 +63,26 @@ TIME_DISCRE_TURB= EULER_IMPLICIT CFL_NUMBER= 10 CFL_REDUCTION_TURB= 1.0 CFL_ADAPT= YES -CFL_ADAPT_PARAM= ( 1, 1.05, 10, 150 ) +CFL_ADAPT_PARAM= ( 1, 1.05, 10, 400 ) % LINEAR_SOLVER= FGMRES LINEAR_SOLVER_PREC= ILU LINEAR_SOLVER_ERROR= 0.2 -LINEAR_SOLVER_ITER= 6 +LINEAR_SOLVER_ITER= 8 % MGLEVEL= 0 NEWTON_KRYLOV= YES NEWTON_KRYLOV_IPARAM= ( 500, 0, 1 ) % n0, np, ft -NEWTON_KRYLOV_DPARAM= ( 0.0, 1e-20, -3, 1e-5, -1 ) % r0, tp, rf, e +NEWTON_KRYLOV_DPARAM= ( 0.0, 1e-20, -3, 1e-5, -1 ) % r0, tp, rf, e, rnk % % ------------------------ CONVERGENCE CRITERIA ------------------------- % % -ITER= 20000 +ITER= 2000 CONV_RESIDUAL_MINVAL= -11.5 +CONV_CAUCHY_ELEMS= 250 +CONV_CAUCHY_EPS= 1e-5 +CONV_FIELD= DRAG, LIFT + % % --------------------------- INPUT / OUTPUT ---------------------------- % % diff --git a/config_template.cfg b/config_template.cfg index bd7a8ff17003..3a31be6f179c 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1645,10 +1645,11 @@ TIME_DISCRE_FLOW= EULER_IMPLICIT NEWTON_KRYLOV= NO % % Integer parameters {startup iters, precond iters, initial tolerance relaxation}. -NEWTON_KRYLOV_IPARAM= (10, 3, 2) +NEWTON_KRYLOV_IPARAM= (20, 3, 2) % -% Double parameters {startup residual drop, precond tolerance, full tolerance residual drop, findiff step}. -NEWTON_KRYLOV_DPARAM= (1.0, 0.1, -6.0, 1e-5) +% Double parameters {startup residual drop, precond tolerance, full tolerance residual drop, findiff step, +% and MUSCL relaxation for the finite difference products}. +NEWTON_KRYLOV_DPARAM= (-2.0, 0.1, -6.0, 1e-5, 1.0) % ------------------- FEM FLOW NUMERICAL METHOD DEFINITION --------------------% % From cee6d239934f72184810ee29ebef9246280cf99b Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 29 Sep 2025 09:08:43 -0700 Subject: [PATCH 4/9] fix merge --- Common/src/CConfig.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 96c319151235..542f79b4b815 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1370,7 +1370,7 @@ void CConfig::SetConfig_Options() { addBoolOption("INC_ENERGY_EQUATION", Energy_Equation, false); /*!\brief TEMPERATURE_LIMITS \n DESCRIPTION: Temperature limits for incompressible flows (0.0, 5000 K by default) \ingroup Config*/ TemperatureLimits[0] = 0.0; TemperatureLimits[1] = 5000.0; - addDoubleArrayOption("TEMPERATURE_LIMITS", 2, TemperatureLimits); + addDoubleArrayOption("TEMPERATURE_LIMITS", 2, false, TemperatureLimits); /*!\brief INC_DENSITY_REF \n DESCRIPTION: Reference density for incompressible flows \ingroup Config*/ addDoubleOption("INC_DENSITY_REF", Inc_Density_Ref, 1.0); /*!\brief INC_VELOCITY_REF \n DESCRIPTION: Reference velocity for incompressible flows (1.0 by default) \ingroup Config*/ From f98e478ed84c529fe005cf7c669d5248b9434d9d Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 29 Sep 2025 19:04:46 -0700 Subject: [PATCH 5/9] fix --- Common/include/option_structure.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/include/option_structure.inl b/Common/include/option_structure.inl index 5dba5c845b22..5f4a6e4fab51 100644 --- a/Common/include/option_structure.inl +++ b/Common/include/option_structure.inl @@ -245,7 +245,7 @@ class COptionArray final : public COptionBase { string SetValue(const vector& option_value) override { COptionBase::SetValue(option_value); // Check that the size is correct - if ((option_value.size() != size_t(size) && allow_fewer) || option_value.size() > size_t(size)) { + if ((option_value.size() < size_t(size) && !allow_fewer) || option_value.size() > size_t(size)) { string newstring; newstring.append(this->name); newstring.append(": wrong number of arguments: "); From 07494fc46055d0ab6727a3838ca6f0e5ee947562 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 4 Oct 2025 12:47:05 -0700 Subject: [PATCH 6/9] reduce number of gradients and limiters for Roe --- .../numerics_simd/flow/convection/common.hpp | 7 +++---- .../numerics_simd/flow/convection/roe.hpp | 5 +++-- SU2_CFD/include/variables/CEulerVariable.hpp | 5 +++++ SU2_CFD/src/solvers/CEulerSolver.cpp | 3 +-- SU2_CFD/src/variables/CEulerVariable.cpp | 16 ++++++++++++++-- 5 files changed, 26 insertions(+), 10 deletions(-) diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index 5b725fd713e1..e681136a3d48 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -114,7 +114,9 @@ FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int const CPair& V1st, const VectorDbl& vector_ij, const VariableType& solution) { - static_assert(ReconVarType::nVar <= PrimVarType::nVar,""); + /*--- Recompute density and enthalpy instead of reconstructing. ---*/ + constexpr auto nVarGrad = ReconVarType::nVar - 2; + static_assert(nVarGrad <= PrimVarType::nVar); const auto& gradients = solution.GetGradient_Reconstruction(); const auto& limiters = solution.GetLimiter_Primitive(); @@ -127,9 +129,6 @@ FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int } if (muscl) { - /*--- Recompute density and enthalpy instead of reconstructing. ---*/ - constexpr auto nVarGrad = ReconVarType::nVar - 2; - switch (limiterType) { case LIMITER::NONE: musclUnlimited(iPoint, vector_ij, 0.5, gradients, V.i.all); diff --git a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp index d46cacea3766..ef53ae64619b 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp @@ -51,7 +51,7 @@ class CRoeBase : public Base { protected: using Base::nDim; static constexpr size_t nVar = CCompressibleConservatives::nVar; - static constexpr size_t nPrimVarGrad = nDim+4; + static constexpr size_t nPrimVarGrad = nDim+2; static constexpr size_t nPrimVar = Max(Base::nPrimVar, nPrimVarGrad); const su2double kappa; @@ -123,7 +123,8 @@ class CRoeBase : public Base { for (size_t iDim = 0; iDim < nDim; ++iDim) { mod_vector_ij(iDim) = nkRelax * vector_ij(iDim); } - auto V = reconstructPrimitives >( + /*--- Recompute density and enthalpy instead of reconstructing. ---*/ + auto V = reconstructPrimitives >( iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, mod_vector_ij, solution); /*--- Compute conservative variables. ---*/ diff --git a/SU2_CFD/include/variables/CEulerVariable.hpp b/SU2_CFD/include/variables/CEulerVariable.hpp index c94e5138590e..4d49f3825721 100644 --- a/SU2_CFD/include/variables/CEulerVariable.hpp +++ b/SU2_CFD/include/variables/CEulerVariable.hpp @@ -30,6 +30,11 @@ #include #include "CFlowVariable.hpp" +/*! + * \brief Returns the number of primitive variables for which to compute gradients. + */ +unsigned long EulerNPrimVarGrad(const CConfig *config, unsigned long ndim); + /*! * \class CEulerVariable * \brief Class for defining the variables of the compressible Euler solver. diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 9e0813f73b19..5b85687071ff 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -64,7 +64,6 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING); const bool adjoint = config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint(); - const bool centered = config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED; int Unst_RestartIter = 0; unsigned long iPoint, iMarker, counter_local = 0, counter_global = 0; @@ -120,7 +119,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, nVar = nDim + 2; nPrimVar = nDim + 9; /*--- Centered schemes only need gradients for viscous fluxes (T and v). ---*/ - nPrimVarGrad = nDim + (centered && !config->GetContinuous_Adjoint() ? 1 : 4); + nPrimVarGrad = EulerNPrimVarGrad(config, nDim); nSecondaryVar = nSecVar; nSecondaryVarGrad = 2; diff --git a/SU2_CFD/src/variables/CEulerVariable.cpp b/SU2_CFD/src/variables/CEulerVariable.cpp index ee29d4992bc6..263530705c87 100644 --- a/SU2_CFD/src/variables/CEulerVariable.cpp +++ b/SU2_CFD/src/variables/CEulerVariable.cpp @@ -28,10 +28,22 @@ #include "../../include/variables/CEulerVariable.hpp" #include "../../include/fluid/CFluidModel.hpp" +unsigned long EulerNPrimVarGrad(const CConfig *config, unsigned long ndim) { + if (config->GetContinuous_Adjoint()) return ndim + 4; + if (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) return ndim + 1; + + const bool ideal_gas = config->GetKind_FluidModel() == STANDARD_AIR || + config->GetKind_FluidModel() == IDEAL_GAS; + if (ideal_gas && config->GetKind_Upwind_Flow() == UPWIND::ROE && !config->Low_Mach_Correction()) { + // Based on CRoeBase (numerics_simd). + return ndim + 2; + } + return ndim + 4; +} + CEulerVariable::CEulerVariable(su2double density, const su2double *velocity, su2double energy, unsigned long npoint, unsigned long ndim, unsigned long nvar, const CConfig *config) - : CFlowVariable(npoint, ndim, nvar, ndim + 9, - ndim + (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED && !config->GetContinuous_Adjoint() ? 1 : 4), config), + : CFlowVariable(npoint, ndim, nvar, ndim + 9, EulerNPrimVarGrad(config, ndim), config), indices(ndim, 0) { const bool dual_time = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || From 0c6cd0ece60956412fb922c9f758ba82e8c2314d Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 4 Oct 2025 21:19:18 -0700 Subject: [PATCH 7/9] fix var counts --- SU2_CFD/include/numerics_simd/flow/convection/common.hpp | 6 +++--- SU2_CFD/include/numerics_simd/flow/convection/roe.hpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index e681136a3d48..dfaf4d5ef8f7 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -114,9 +114,7 @@ FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int const CPair& V1st, const VectorDbl& vector_ij, const VariableType& solution) { - /*--- Recompute density and enthalpy instead of reconstructing. ---*/ - constexpr auto nVarGrad = ReconVarType::nVar - 2; - static_assert(nVarGrad <= PrimVarType::nVar); + static_assert(ReconVarType::nVar <= PrimVarType::nVar); const auto& gradients = solution.GetGradient_Reconstruction(); const auto& limiters = solution.GetLimiter_Primitive(); @@ -129,6 +127,8 @@ FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int } if (muscl) { + /*--- Recompute density and enthalpy instead of reconstructing. ---*/ + constexpr auto nVarGrad = ReconVarType::nVar - 2; switch (limiterType) { case LIMITER::NONE: musclUnlimited(iPoint, vector_ij, 0.5, gradients, V.i.all); diff --git a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp index ef53ae64619b..2460267aaa87 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/roe.hpp @@ -51,7 +51,7 @@ class CRoeBase : public Base { protected: using Base::nDim; static constexpr size_t nVar = CCompressibleConservatives::nVar; - static constexpr size_t nPrimVarGrad = nDim+2; + static constexpr size_t nPrimVarGrad = nDim+4; static constexpr size_t nPrimVar = Max(Base::nPrimVar, nPrimVarGrad); const su2double kappa; @@ -124,7 +124,7 @@ class CRoeBase : public Base { mod_vector_ij(iDim) = nkRelax * vector_ij(iDim); } /*--- Recompute density and enthalpy instead of reconstructing. ---*/ - auto V = reconstructPrimitives >( + auto V = reconstructPrimitives >( iEdge, iPoint, jPoint, gamma, gasConst, muscl, typeLimiter, V1st, mod_vector_ij, solution); /*--- Compute conservative variables. ---*/ From f2c50ccaee15c76424f66ac4ce843bb9ecdf3db9 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 5 Oct 2025 08:27:32 -0700 Subject: [PATCH 8/9] update a few tests --- SU2_CFD/src/output/CFlowCompOutput.cpp | 8 ++++++-- TestCases/parallel_regression.py | 2 +- TestCases/serial_regression_AD.py | 2 +- TestCases/vandv.py | 5 ++--- 4 files changed, 10 insertions(+), 7 deletions(-) diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index a765ea72c43a..2476498e84d7 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -277,6 +277,7 @@ void CFlowCompOutput::SetVolumeOutputFields(CConfig *config){ SetVolumeOutputFieldsScalarResidual(config); if (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE && config->GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE) { + AddVolumeOutput("LIMITER_TEMPERATURE", "Limiter_Temperature", "LIMITER", "Limiter value of the temperature"); AddVolumeOutput("LIMITER_VELOCITY-X", "Limiter_Velocity_x", "LIMITER", "Limiter value of the x-velocity"); AddVolumeOutput("LIMITER_VELOCITY-Y", "Limiter_Velocity_y", "LIMITER", "Limiter value of the y-velocity"); if (nDim == 3) { @@ -364,14 +365,17 @@ void CFlowCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv } if (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE && config->GetKind_SlopeLimit_Flow() != LIMITER::VAN_ALBADA_EDGE) { + SetVolumeOutputValue("LIMITER_TEMPERATURE", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, 0)); SetVolumeOutputValue("LIMITER_VELOCITY-X", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, 1)); SetVolumeOutputValue("LIMITER_VELOCITY-Y", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, 2)); if (nDim == 3){ SetVolumeOutputValue("LIMITER_VELOCITY-Z", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, 3)); } SetVolumeOutputValue("LIMITER_PRESSURE", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, nDim+1)); - SetVolumeOutputValue("LIMITER_DENSITY", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, nDim+2)); - SetVolumeOutputValue("LIMITER_ENTHALPY", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, nDim+3)); + if (solver[FLOW_SOL]->GetnPrimVarGrad() > nDim + 2) { + SetVolumeOutputValue("LIMITER_DENSITY", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, nDim+2)); + SetVolumeOutputValue("LIMITER_ENTHALPY", iPoint, Node_Flow->GetLimiter_Primitive(iPoint, nDim+3)); + } } if (config->GetKind_RoeLowDiss() != NO_ROELOWDISS){ diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 86a534b6ea1f..b32e24515364 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -431,7 +431,7 @@ def main(): turb_oneram6_nk.cfg_dir = "rans/oneram6" turb_oneram6_nk.cfg_file = "turb_ONERAM6_nk.cfg" turb_oneram6_nk.test_iter = 20 - turb_oneram6_nk.test_vals = [-4.827571, -4.425650, -11.379658, 0.224787, 0.044208, 1.000000, -0.642711, 31.384000] + turb_oneram6_nk.test_vals = [-5.262975, -4.885414, -11.509429, 0.218369, 0.067725, 2, -0.772645, 10] turb_oneram6_nk.timeout = 600 turb_oneram6_nk.tol = 0.0001 test_list.append(turb_oneram6_nk) diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 268a1cd78c4a..7fc8224b2fca 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -205,7 +205,7 @@ def main(): discadj_fsi.cfg_dir = "disc_adj_fsi" discadj_fsi.cfg_file = "config.cfg" discadj_fsi.test_iter = 6 - discadj_fsi.test_vals = [6.000000, -8.929426, -10.024862, 0.000000, -0.000002] + discadj_fsi.test_vals = [6, -8.927769, -10.148808, 3.1045e-11, -1.7610e-06] test_list.append(discadj_fsi) ################################### diff --git a/TestCases/vandv.py b/TestCases/vandv.py index 20aab18891c6..3aedecc2b3c9 100644 --- a/TestCases/vandv.py +++ b/TestCases/vandv.py @@ -45,8 +45,7 @@ def main(): p30n30.cfg_dir = "vandv/rans/30p30n" p30n30.cfg_file = "config.cfg" p30n30.test_iter = 5 - p30n30.test_vals = [-11.498923, -11.508761, -11.978899, -11.700024, -14.235389, 0.052235, 2.830394, 1.318894, -0.289257] - p30n30.test_vals_aarch64 = [-11.498923, -11.508761, -11.978899, -11.700024, -14.235389, 0.052235, 2.830394, 1.318894, -0.289257] + p30n30.test_vals = [-11.267106, -11.168215, -11.182822, -10.949673, -14.233489, 0.052235, 2.830394, 1.318894, -1.210648, 1, 1.2763e+01] test_list.append(p30n30) # flat plate - sst-v1994m @@ -85,7 +84,7 @@ def main(): swbli_sst.test_vals = [-12.001545, -11.350636, -12.056760, -10.870102, -11.411568, -2.263450, 0.001796, -1.450519, -2.930524, 10.000000] test_list.append(swbli_sst) - # DSMA661 - SA + # DSMA661 - SA dsma661_sa = TestCase('dsma661_sa') dsma661_sa.cfg_dir = "vandv/rans/dsma661" dsma661_sa.cfg_file = "dsma661_sa_config.cfg" From a5f5513a38a009b076538382b90a58704072e276 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 5 Oct 2025 15:45:55 -0700 Subject: [PATCH 9/9] same logic for target tol --- SU2_CFD/src/integration/CNewtonIntegration.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/integration/CNewtonIntegration.cpp b/SU2_CFD/src/integration/CNewtonIntegration.cpp index a3286fd7c012..e7f0416e623c 100644 --- a/SU2_CFD/src/integration/CNewtonIntegration.cpp +++ b/SU2_CFD/src/integration/CNewtonIntegration.cpp @@ -277,9 +277,10 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver ** solvers[FLOW_SOL]->SetResLinSolver(eps); if (!startupPeriod && autoRelaxation) { - if (eps > config->GetLinear_Solver_Error()) { + const su2double adaptTol = config->GetCFL_Adapt() ? config->GetCFL_AdaptParam(4) : 0; + if (eps > fmax(config->GetLinear_Solver_Error(), adaptTol)) { nkRelaxation *= 0.9; - } else if (iter < config->GetLinear_Solver_Iter()) { + } else if (eps < 0.9 * fmax(config->GetLinear_Solver_Error(), adaptTol)) { nkRelaxation = fmin(nkRelaxation * 1.05, 1); } }