diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index abbb0f401da6..c996c643e1f1 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); @@ -4375,12 +4376,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..5f4a6e4fab51 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 a2a529917bf3..542f79b4b815 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); @@ -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*/ @@ -1381,7 +1381,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*/ @@ -1406,10 +1406,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); @@ -1429,7 +1429,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: */ @@ -1626,7 +1626,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); @@ -1640,13 +1640,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*/ @@ -1654,7 +1654,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*/ @@ -1727,14 +1727,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); @@ -1756,7 +1756,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); @@ -1779,9 +1779,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()); - /* DESCRIPTION: Double parameters {startup residual drop, precond tolerance, full tolerance residual drop, findiff step}. */ - addDoubleArrayOption("NEWTON_KRYLOV_DPARAM", NK_DblParam.size(), NK_DblParam.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, NK relaxation}. */ + 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); @@ -1989,7 +1989,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*/ @@ -2010,7 +2010,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); @@ -2076,7 +2076,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 */ @@ -2088,7 +2088,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) */ @@ -2166,11 +2166,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); @@ -2260,21 +2260,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) */ @@ -2356,7 +2356,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); @@ -2405,7 +2405,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 ---*/ @@ -2524,7 +2524,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); @@ -2570,7 +2570,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); @@ -2709,10 +2709,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 */ @@ -2772,7 +2772,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); @@ -2812,7 +2812,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/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index 5b725fd713e1..dfaf4d5ef8f7 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -114,7 +114,7 @@ FORCEINLINE CPair reconstructPrimitives(Int iEdge, Int iPoint, Int const CPair& V1st, const VectorDbl& vector_ij, const VariableType& solution) { - static_assert(ReconVarType::nVar <= PrimVarType::nVar,""); + static_assert(ReconVarType::nVar <= PrimVarType::nVar); const auto& gradients = solution.GetGradient_Reconstruction(); const auto& limiters = solution.GetLimiter_Primitive(); @@ -129,7 +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 7309c0f3f678..2460267aaa87 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,13 @@ 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); + } + /*--- Recompute density and enthalpy instead of reconstructing. ---*/ 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/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/integration/CNewtonIntegration.cpp b/SU2_CFD/src/integration/CNewtonIntegration.cpp index daaa7fcd64fc..e7f0416e623c 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 && startupIters > 0 && !config->GetRestart()) { + 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,15 @@ 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) { + const su2double adaptTol = config->GetCFL_Adapt() ? config->GetCFL_AdaptParam(4) : 0; + if (eps > fmax(config->GetLinear_Solver_Error(), adaptTol)) { + nkRelaxation *= 0.9; + } else if (eps < 0.9 * fmax(config->GetLinear_Solver_Error(), adaptTol)) { + nkRelaxation = fmin(nkRelaxation * 1.05, 1); + } + } } END_SU2_OMP_SAFE_GLOBAL_ACCESS @@ -303,11 +340,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/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/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 73dc2a6189c5..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; @@ -1796,6 +1795,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 +1875,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/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) || 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/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/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/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" diff --git a/TestCases/vandv/rans/30p30n/config.cfg b/TestCases/vandv/rans/30p30n/config.cfg index da6749b88b36..f4b9fa384a37 100644 --- a/TestCases/vandv/rans/30p30n/config.cfg +++ b/TestCases/vandv/rans/30p30n/config.cfg @@ -60,24 +60,29 @@ 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, 400 ) % LINEAR_SOLVER= FGMRES LINEAR_SOLVER_PREC= ILU LINEAR_SOLVER_ERROR= 0.2 -LINEAR_SOLVER_ITER= 4 +LINEAR_SOLVER_ITER= 8 % 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, rnk % % ------------------------ CONVERGENCE CRITERIA ------------------------- % % -ITER= 20 +ITER= 2000 CONV_RESIDUAL_MINVAL= -11.5 +CONV_CAUCHY_ELEMS= 250 +CONV_CAUCHY_EPS= 1e-5 +CONV_FIELD= DRAG, LIFT + % % --------------------------- INPUT / OUTPUT ---------------------------- % % @@ -89,5 +94,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 ) diff --git a/config_template.cfg b/config_template.cfg index 2b43a0646784..ef0d921c3ffd 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -1648,10 +1648,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 --------------------% %