diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 2f900df73aae..626b8c30f9a3 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -150,7 +150,6 @@ class CConfig { su2double CL_Target; /*!< \brief Fixed Cl mode Target Cl. */ su2double Confinement_Param; /*!< \brief Confinement paramenter for Vorticity Confinement method. */ TIME_MARCHING TimeMarching; /*!< \brief Steady or unsteady (time stepping or dual time stepping) computation. */ - unsigned short Dynamic_Analysis; /*!< \brief Static or dynamic structural analysis. */ su2double FixAzimuthalLine; /*!< \brief Fix an azimuthal line due to misalignments of the nearfield. */ su2double **DV_Value; /*!< \brief Previous value of the design variable. */ su2double Venkat_LimiterCoeff; /*!< \brief Limiter coefficient */ @@ -168,9 +167,6 @@ class CConfig { su2double HarmonicBalance_Period; /*!< \brief Period of oscillation to be used with harmonic balance computations. */ su2double Delta_UnstTime, /*!< \brief Time step for unsteady computations. */ Delta_UnstTimeND; /*!< \brief Time step for unsteady computations (non dimensional). */ - su2double Delta_DynTime, /*!< \brief Time step for dynamic structural computations. */ - Total_DynTime, /*!< \brief Total time for dynamic structural computations. */ - Current_DynTime; /*!< \brief Global time of the dynamic structural computations. */ su2double Total_UnstTime, /*!< \brief Total time for unsteady computations. */ Total_UnstTimeND; /*!< \brief Total time for unsteady computations (non dimensional). */ su2double Current_UnstTime, /*!< \brief Global time of the unsteady simulation. */ @@ -2112,17 +2108,11 @@ class CConfig { su2double GetElasticyMod(unsigned short id_val) const { return ElasticityMod[id_val]; } /*! - * \brief Decide whether to apply DE effects to the model. - * \return TRUE if the DE effects are to be applied, FALSE otherwise. - */ + * \brief Decide whether to apply DE effects to the model. + * \return TRUE if the DE effects are to be applied, FALSE otherwise. + */ bool GetDE_Effects(void) const { return DE_Effects; } - /*! - * \brief Decide whether to predict the DE effects for the next time step. - * \return TRUE if the DE effects are to be applied, FALSE otherwise. - */ - bool GetDE_Predicted(void); - /*! * \brief Get the number of different electric constants. * \return Value of the DE modulus. @@ -8722,34 +8712,6 @@ class CConfig { */ bool GetSteadyRestart(void) const { return SteadyRestart; } - /*! - * \brief Provides information about the time integration of the structural analysis, and change the write in the output - * files information about the iteration. - * \return The kind of time integration: Static or dynamic analysis - */ - unsigned short GetDynamic_Analysis(void) const { return Dynamic_Analysis; } - - /*! - * \brief If we are prforming an unsteady simulation, there is only - * one value of the time step for the complete simulation. - * \return Value of the time step in an unsteady simulation (non dimensional). - */ - su2double GetDelta_DynTime(void) const { return Delta_DynTime; } - - /*! - * \brief If we are prforming an unsteady simulation, there is only - * one value of the time step for the complete simulation. - * \return Value of the time step in an unsteady simulation (non dimensional). - */ - su2double GetTotal_DynTime(void) const { return Total_DynTime; } - - /*! - * \brief If we are prforming an unsteady simulation, there is only - * one value of the time step for the complete simulation. - * \return Value of the time step in an unsteady simulation (non dimensional). - */ - su2double GetCurrent_DynTime(void) const { return Current_DynTime; } - /*! * \brief Get the current instance. * \return Current instance identifier. @@ -8781,22 +8743,10 @@ class CConfig { unsigned short GetnIntCoeffs(void) const { return nIntCoeffs; } /*! - * \brief Get the number of different values for the elasticity modulus. - * \return Number of different values for the elasticity modulus. + * \brief Get the number of different materials for the elasticity solver. + * \return Number of different materials. */ - unsigned short GetnElasticityMod(void) const { return nElasticityMod; } - - /*! - * \brief Get the number of different values for the Poisson ratio. - * \return Number of different values for the Poisson ratio. - */ - unsigned short GetnPoissonRatio(void) const { return nPoissonRatio; } - - /*! - * \brief Get the number of different values for the Material density. - * \return Number of different values for the Material density. - */ - unsigned short GetnMaterialDensity(void) const { return nMaterialDensity; } + unsigned short GetnElasticityMat(void) const { return nElasticityMod; } /*! * \brief Get the integration coefficients for the Generalized Alpha - Newmark integration integration scheme. @@ -9282,12 +9232,6 @@ class CConfig { */ void SetnTime_Iter(unsigned long val_iter) { nTimeIter = val_iter; } - /*! - * \brief Get the number of pseudo-time iterations - * \return Number of pseudo-time steps run for the single-zone problem - */ - unsigned long GetnIter(void) const { return nIter; } - /*! * \brief Get the restart iteration * \return Iteration for the restart of multizone problems diff --git a/Common/include/geometry/elements/CElementProperty.hpp b/Common/include/geometry/elements/CElementProperty.hpp index fce7425694e0..90efe7194b72 100644 --- a/Common/include/geometry/elements/CElementProperty.hpp +++ b/Common/include/geometry/elements/CElementProperty.hpp @@ -93,7 +93,7 @@ class CProperty { /*! * \brief Extract the derivative of the Design density. */ - inline virtual su2double GetAdjointDensity(void) const { return 0.0; } + inline virtual su2double GetAdjointDensity(void) { return 0.0; } /*! * \brief Register the Design density as an AD input variable. @@ -177,7 +177,11 @@ class CElementProperty final : public CProperty { /*! * \brief Extract the derivative of the Design density. */ - inline su2double GetAdjointDensity(void) const override { return SU2_TYPE::GetDerivative(design_rho); } + inline su2double GetAdjointDensity(void) override { + su2double der = SU2_TYPE::GetDerivative(design_rho); + AD::ResetInput(design_rho); + return der; + } /*! * \brief Register the Design density as an AD input variable. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 954088472a7e..69acd1c0c671 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1434,14 +1434,12 @@ static const MapType Heat_TimeStep_Map = { * \brief Type of time integration schemes */ enum class STRUCT_TIME_INT { - CD_EXPLICIT, /*!< \brief Support for implementing an explicit method. */ NEWMARK_IMPLICIT, /*!< \brief Implicit Newmark integration definition. */ GENERALIZED_ALPHA, /*!< \brief Support for implementing another implicit method. */ }; static const MapType Time_Int_Map_FEA = { - MakePair("CD_EXPLICIT", STRUCT_TIME_INT::CD_EXPLICIT) MakePair("NEWMARK_IMPLICIT", STRUCT_TIME_INT::NEWMARK_IMPLICIT) - MakePair("GENERALIZED_ALPHA", STRUCT_TIME_INT::GENERALIZED_ALPHA) + // MakePair("GENERALIZED_ALPHA", STRUCT_TIME_INT::GENERALIZED_ALPHA) Not fully implemented. }; /*! @@ -2368,18 +2366,6 @@ enum class RECORDING { SOLUTION_AND_MESH, }; -/*! - * \brief Types of schemes for dynamic structural computations - */ -enum ENUM_DYNAMIC { - STATIC = 0, /*!< \brief A static structural computation. */ - DYNAMIC = 1 /*!< \brief Use a time stepping strategy for dynamic computations. */ -}; -static const MapType Dynamic_Map = { - MakePair("NO", STATIC) - MakePair("YES", DYNAMIC) -}; - /*! * \brief Types of input file formats */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 54ffa66ae07f..e705b7f38d40 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2435,12 +2435,6 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Temporary: pseudo static analysis (no density in dynamic analysis) * Options: NO, YES \ingroup Config */ addBoolOption("PSEUDO_STATIC", PseudoStatic, false); - /* DESCRIPTION: Dynamic or static structural analysis */ - addEnumOption("DYNAMIC_ANALYSIS", Dynamic_Analysis, Dynamic_Map, STATIC); - /* DESCRIPTION: Time Step for dynamic analysis (s) */ - addDoubleOption("DYN_TIMESTEP", Delta_DynTime, 0.0); - /* DESCRIPTION: Total Physical Time for dual time stepping simulations (s) */ - addDoubleOption("DYN_TIME", Total_DynTime, 1.0); /* DESCRIPTION: Parameter alpha for Newmark scheme (s) */ addDoubleOption("NEWMARK_BETA", Newmark_beta, 0.25); /* DESCRIPTION: Parameter delta for Newmark scheme (s) */ @@ -2987,6 +2981,12 @@ void CConfig::SetConfig_Parsing(istream& config_buffer){ newString.append("\n"); if (!option_name.compare("SINGLEZONE_DRIVER")) newString.append("Option SINGLEZONE_DRIVER is deprecated, it does not have a replacement.\n\n"); + else if (!option_name.compare("DYN_TIMESTEP")) + newString.append("DYN_TIMESTEP is deprecated. Use TIME_STEP instead.\n\n"); + else if (!option_name.compare("DYN_TIME")) + newString.append("DYN_TIME is deprecated. Use MAX_TIME instead.\n\n"); + else if (!option_name.compare("DYNAMIC_ANALYSIS")) + newString.append("DYNAMIC_ANALYSIS is deprecated. Use TIME_DOMAIN instead.\n\n"); else { /*--- Find the most likely candidate for the unrecognized option, based on the length of start and end character sequences shared by candidates and the option. ---*/ @@ -3660,7 +3660,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i if (Time_Domain){ Delta_UnstTime = Time_Step; - Delta_DynTime = Time_Step; if (TimeMarching == TIME_MARCHING::TIME_STEPPING){ InnerIter = 1; } @@ -4649,6 +4648,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i MaterialDensity = new su2double[1]; MaterialDensity[0] = 7854; } + if (nElasticityMod != nPoissonRatio || nElasticityMod != nMaterialDensity) { + SU2_MPI::Error("ELASTICITY_MODULUS, POISSON_RATIO, and MATERIAL_DENSITY need to have the same number " + "of entries (the number of materials).", CURRENT_FUNCTION); + } + if (nElectric_Constant == 0) { nElectric_Constant = 1; Electric_Constant = new su2double[1]; Electric_Constant[0] = 0.0; @@ -6803,7 +6807,7 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { else { if (Time_Domain) { cout << "Dynamic structural analysis."<< endl; - cout << "Time step provided by the user for the dynamic analysis(s): "<< Delta_DynTime << "." << endl; + cout << "Time step provided by the user for the dynamic analysis(s): "<< Time_Step << "." << endl; } else { cout << "Static structural analysis." << endl; } @@ -6871,14 +6875,11 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { if (fea) { switch (Kind_TimeIntScheme_FEA) { - case STRUCT_TIME_INT::CD_EXPLICIT: - cout << "Explicit time integration (NOT IMPLEMENTED YET)." << endl; - break; case STRUCT_TIME_INT::GENERALIZED_ALPHA: cout << "Generalized-alpha method." << endl; break; case STRUCT_TIME_INT::NEWMARK_IMPLICIT: - if (Dynamic_Analysis) cout << "Newmark implicit method for the structural time integration." << endl; + if (Time_Domain) cout << "Newmark implicit method for the structural time integration." << endl; switch (Kind_Linear_Solver) { case BCGSTAB: cout << "BCGSTAB is used for solving the linear system." << endl; @@ -8416,9 +8417,6 @@ void CConfig::SetGlobalParam(MAIN_SOLVER val_solver, case MAIN_SOLVER::FEM_ELASTICITY: case MAIN_SOLVER::DISC_ADJ_FEM: - - Current_DynTime = static_cast(TimeIter)*Delta_DynTime; - if (val_system == RUNTIME_FEA_SYS) { SetKind_ConvNumScheme(NONE, CENTERED::NONE, UPWIND::NONE, LIMITER::NONE, NONE, NONE); SetKind_TimeIntScheme(NONE); @@ -9994,7 +9992,6 @@ void CConfig::SetMultizone(const CConfig *driver_config, const CConfig* const* c /*--- Fix the Time Step for all subdomains, for the case of time-dependent problems ---*/ if (driver_config->GetTime_Domain()){ Delta_UnstTime = driver_config->GetTime_Step(); - Delta_DynTime = driver_config->GetTime_Step(); Time_Domain = true; } diff --git a/SU2_CFD/include/drivers/CDriverBase.hpp b/SU2_CFD/include/drivers/CDriverBase.hpp index 4e79537a8298..d720246efb74 100644 --- a/SU2_CFD/include/drivers/CDriverBase.hpp +++ b/SU2_CFD/include/drivers/CDriverBase.hpp @@ -499,6 +499,15 @@ class CDriverBase { main_geometry->GetnVertex(iMarker), "MarkerPrimitives", false); } + /*! + * \brief Get a read-only view of the geometry sensitivity of a discrete adjoint solver. + */ + inline CPyWrapperMatrixView Sensitivity(unsigned short iSolver) { + auto* solver = GetSolverAndCheckMarker(iSolver); + auto& sensitivity = const_cast(solver->GetNodes()->GetSensitivity()); + return CPyWrapperMatrixView(sensitivity, "Sensitivity", true); + } + /*! * \brief Set the temperature of a vertex on a specified marker (MARKER_PYTHON_CUSTOM). * \note This can be the input of a heat or flow solver in a CHT setting. diff --git a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp index abd8fcddfb0d..2755731f5a7b 100644 --- a/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp +++ b/SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp @@ -50,7 +50,8 @@ class CDiscAdjFEASolver final : public CSolver { su2double* val = nullptr; /*!< \brief Value of the variable. */ su2double* LocalSens = nullptr; /*!< \brief Local sensitivity (domain). */ su2double* GlobalSens = nullptr; /*!< \brief Global sensitivity (mpi). */ - su2double* TotalSens = nullptr; /*!< \brief Total sensitivity (time domain). */ + su2double* OldSens = nullptr; /*!< \brief Previous global sensitivity, used to update the total. */ + su2double* TotalSens = nullptr; /*!< \brief Total sensitivity (integrated over time). */ su2double& operator[] (unsigned short i) { return val[i]; } const su2double& operator[] (unsigned short i) const { return val[i]; } @@ -61,6 +62,7 @@ class CDiscAdjFEASolver final : public CSolver { val = new su2double[n](); LocalSens = new su2double[n](); GlobalSens = new su2double[n](); + OldSens = new su2double[n](); TotalSens = new su2double[n](); } @@ -69,6 +71,7 @@ class CDiscAdjFEASolver final : public CSolver { delete [] val; delete [] LocalSens; delete [] GlobalSens; + delete [] OldSens; delete [] TotalSens; } @@ -80,10 +83,19 @@ class CDiscAdjFEASolver final : public CSolver { for (auto i = 0u; i < size; ++i) LocalSens[i] = SU2_TYPE::GetDerivative(val[i]); SU2_MPI::Allreduce(LocalSens, GlobalSens, size, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + + for (auto i = 0u; i < size; ++i) { + /*--- Update the total by subtracting the old and adding the new value. + * Then update the old value for the next call to this function. ---*/ + TotalSens[i] += GlobalSens[i] - OldSens[i]; + OldSens[i] = GlobalSens[i]; + } } - void UpdateTotal() { - for (auto i = 0u; i < size; ++i) TotalSens[i] += GlobalSens[i]; + void Store() { + /*--- Clears the old values such that on the next time step the total is + * incremented instead of updated. ---*/ + for (auto i = 0u; i < size; ++i) OldSens[i] = 0.0; } ~SensData() { clear(); } @@ -213,42 +225,6 @@ class CDiscAdjFEASolver final : public CSolver { */ inline su2double GetTotal_Sens_DVFEA(unsigned short iDVFEA) const override { return DV.TotalSens[iDVFEA]; } - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the Young Modulus E - */ - inline su2double GetGlobal_Sens_E(unsigned short iVal) const override { return E.GlobalSens[iVal]; } - - /*! - * \brief A virtual member. - * \return Value of the Mach sensitivity for the Poisson's ratio Nu - */ - inline su2double GetGlobal_Sens_Nu(unsigned short iVal) const override { return Nu.GlobalSens[iVal]; } - - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the Electric Field in the region iEField - */ - inline su2double GetGlobal_Sens_EField(unsigned short iEField) const override { return EField.GlobalSens[iEField]; } - - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the FEA DV in the region iDVFEA - */ - inline su2double GetGlobal_Sens_DVFEA(unsigned short iDVFEA) const override { return DV.GlobalSens[iDVFEA]; } - - /*! - * \brief Get the total sensitivity for the structural density - * \return Value of the structural density sensitivity - */ - inline su2double GetGlobal_Sens_Rho(unsigned short iVal) const override { return Rho.GlobalSens[iVal]; } - - /*! - * \brief Get the total sensitivity for the structural weight - * \return Value of the structural weight sensitivity - */ - inline su2double GetGlobal_Sens_Rho_DL(unsigned short iVal) const override { return Rho_DL.GlobalSens[iVal]; } - /*! * \brief Get the value of the Young modulus from the adjoint solver * \return Value of the Young modulus from the adjoint solver diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index ed0aee60236d..86eb0813e6a5 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -682,28 +682,9 @@ class CFEASolver : public CFEASolverBase { inline su2double GetFSI_ConvValue(unsigned short val_index) const final { return FSI_Conv[val_index]; } /*! - * \brief Retrieve the value of the dynamic Aitken relaxation factor. - * \return Value of the dynamic Aitken relaxation factor. + * \brief Store the value of the last Aitken relaxation factor in the current time step. */ - inline su2double GetWAitken_Dyn(void) const final { return WAitken_Dyn; } - - /*! - * \brief Retrieve the value of the last Aitken relaxation factor in the previous time step. - * \return Value of the last Aitken relaxation factor in the previous time step. - */ - inline su2double GetWAitken_Dyn_tn1(void) const final { return WAitken_Dyn_tn1; } - - /*! - * \brief Set the value of the dynamic Aitken relaxation factor - * \param[in] Value of the dynamic Aitken relaxation factor - */ - inline void SetWAitken_Dyn(su2double waitk) final { WAitken_Dyn = waitk; } - - /*! - * \brief Set the value of the last Aitken relaxation factor in the current time step. - * \param[in] Value of the last Aitken relaxation factor in the current time step. - */ - inline void SetWAitken_Dyn_tn1(su2double waitk_tn1) final { WAitken_Dyn_tn1 = waitk_tn1; } + inline void SetWAitken_Dyn_tn1() final { WAitken_Dyn_tn1 = WAitken_Dyn; } /*! * \brief Set the value of the load increment for nonlinear structural analysis diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index cdfa19461513..ba28bf8b1933 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3201,42 +3201,6 @@ class CSolver { */ inline virtual su2double GetTotal_Sens_DVFEA(unsigned short iDVFEA) const { return 0.0; } - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the Young Modulus E - */ - inline virtual su2double GetGlobal_Sens_E(unsigned short iVal) const { return 0.0; } - - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the Poisson's ratio Nu - */ - inline virtual su2double GetGlobal_Sens_Nu(unsigned short iVal) const { return 0.0; } - - /*! - * \brief A virtual member. - * \return Value of the structural density sensitivity - */ - inline virtual su2double GetGlobal_Sens_Rho(unsigned short iVal) const { return 0.0; } - - /*! - * \brief A virtual member. - * \return Value of the structural weight sensitivity - */ - inline virtual su2double GetGlobal_Sens_Rho_DL(unsigned short iVal) const { return 0.0; } - - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the Electric Field in the region iEField - */ - inline virtual su2double GetGlobal_Sens_EField(unsigned short iEField) const { return 0.0; } - - /*! - * \brief A virtual member. - * \return Value of the sensitivity coefficient for the FEA DV in the region iDVFEA - */ - inline virtual su2double GetGlobal_Sens_DVFEA(unsigned short iDVFEA) const { return 0.0; } - /*! * \brief A virtual member. * \return Value of the Young modulus from the adjoint solver @@ -3698,27 +3662,8 @@ class CSolver { /*! * \brief A virtual member. - * \return Value of the dynamic Aitken relaxation factor - */ - inline virtual su2double GetWAitken_Dyn(void) const { return 0; } - - /*! - * \brief A virtual member. - * \return Value of the last Aitken relaxation factor in the previous time step. - */ - inline virtual su2double GetWAitken_Dyn_tn1(void) const { return 0; } - - /*! - * \brief A virtual member. - * \param[in] Value of the dynamic Aitken relaxation factor - */ - inline virtual void SetWAitken_Dyn(su2double waitk) { } - - /*! - * \brief A virtual member. - * \param[in] Value of the last Aitken relaxation factor in the previous time step. */ - inline virtual void SetWAitken_Dyn_tn1(su2double waitk_tn1) { } + inline virtual void SetWAitken_Dyn_tn1() { } /*! * \brief A virtual member. diff --git a/SU2_CFD/include/variables/CDiscAdjVariable.hpp b/SU2_CFD/include/variables/CDiscAdjVariable.hpp index 72ca67f9c24a..bc72508de78c 100644 --- a/SU2_CFD/include/variables/CDiscAdjVariable.hpp +++ b/SU2_CFD/include/variables/CDiscAdjVariable.hpp @@ -72,6 +72,7 @@ class CDiscAdjVariable : public CVariable { inline su2double GetSensitivity(unsigned long iPoint, unsigned long iDim) const final { return Sensitivity(iPoint,iDim); } + inline const MatrixType& GetSensitivity() const final { return Sensitivity; } /*! * \brief Set/store the dual time contributions to the adjoint variable. diff --git a/SU2_CFD/include/variables/CFEABoundVariable.hpp b/SU2_CFD/include/variables/CFEABoundVariable.hpp index e04d40467c1b..e0f8c273cc14 100644 --- a/SU2_CFD/include/variables/CFEABoundVariable.hpp +++ b/SU2_CFD/include/variables/CFEABoundVariable.hpp @@ -157,12 +157,12 @@ class CFEABoundVariable final : public CFEAVariable { /*! * \brief Register the flow tractions as input variable. */ - void RegisterFlowTraction() override; + void RegisterFlowTraction(bool reset) override; /*! * \brief Extract the flow traction derivatives. */ - inline su2double ExtractFlowTraction_Sensitivity(unsigned long iPoint, unsigned long iDim) const override { + inline su2double ExtractFlowTractionSensitivity(unsigned long iPoint, unsigned long iDim) const override { if (!fsi_analysis) return 0.0; if (!VertexMap.GetVertexIndex(iPoint)) return 0.0; return SU2_TYPE::GetDerivative(FlowTraction(iPoint,iDim)); diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index da68a28f051d..aa75d07596db 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2160,12 +2160,12 @@ class CVariable { /*! * \brief A virtual member. */ - inline virtual void RegisterFlowTraction() { } + inline virtual void RegisterFlowTraction(bool reset) { } /*! * \brief A virtual member. */ - inline virtual su2double ExtractFlowTraction_Sensitivity(unsigned long iPoint, unsigned long iDim) const { return 0.0; } + inline virtual su2double ExtractFlowTractionSensitivity(unsigned long iPoint, unsigned long iDim) const { return 0.0; } /*! * \brief Register the variables in the solution array as input/output variable. @@ -2224,6 +2224,7 @@ class CVariable { * \return value of the Sensitivity */ inline virtual su2double GetSensitivity(unsigned long iPoint, unsigned long iDim) const { return 0.0; } + inline virtual const MatrixType& GetSensitivity() const { AssertOverride(); return Solution; } inline virtual void SetTau_Wall(unsigned long iPoint, su2double tau_wall) {} diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 586c9cabfdca..6277f409b8d8 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -763,15 +763,15 @@ void CDiscAdjMultizoneDriver::SetObjFunction(RECORDING kind_recording) { } void CDiscAdjMultizoneDriver::SetAdjObjFunction() { - - const auto IterAvg_Obj = config_container[ZONE_0]->GetIter_Avg_Objective(); su2double seeding = 1.0; - if (config_container[ZONE_0]->GetTime_Marching() != TIME_MARCHING::STEADY){ - if (TimeIter < IterAvg_Obj){ - // Default behavior (in case no specific window is chosen) is to use Square-Windowing, i.e. the numerator equals 1.0 + if (config_container[ZONE_0]->GetTime_Domain()) { + const auto IterAvg_Obj = config_container[ZONE_0]->GetIter_Avg_Objective(); + if (TimeIter < IterAvg_Obj) { + /*--- Default behavior when no window is chosen is to use Square-Windowing, i.e. the numerator equals 1.0 ---*/ auto windowEvaluator = CWindowingTools(); - su2double weight = windowEvaluator.GetWndWeight(config_container[ZONE_0]->GetKindWindow(), TimeIter, IterAvg_Obj-1); + const su2double weight = + windowEvaluator.GetWndWeight(config_container[ZONE_0]->GetKindWindow(), TimeIter, IterAvg_Obj - 1); seeding = weight / IterAvg_Obj; } else { diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp index 5a42be040d71..f9cfb269a8ed 100644 --- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp @@ -319,23 +319,21 @@ void CDiscAdjSinglezoneDriver::SetRecording(RECORDING kind_recording){ } void CDiscAdjSinglezoneDriver::SetAdjObjFunction(){ - - const auto IterAvg_Obj = config->GetIter_Avg_Objective(); su2double seeding = 1.0; - CWindowingTools windowEvaluator = CWindowingTools(); - - if (config->GetTime_Marching() != TIME_MARCHING::STEADY){ - if (TimeIter < IterAvg_Obj){ - /*--- Default behavior (in case no specific window is chosen) is to use Square-Windowing, i.e. the numerator equals 1.0 ---*/ - seeding = windowEvaluator.GetWndWeight(config->GetKindWindow(),TimeIter, IterAvg_Obj-1)/ (static_cast(IterAvg_Obj)); + if (config->GetTime_Domain()) { + const auto IterAvg_Obj = config->GetIter_Avg_Objective(); + if (TimeIter < IterAvg_Obj) { + /*--- Default behavior when no window is chosen is to use Square-Windowing, i.e. the numerator equals 1.0 ---*/ + auto windowEvaluator = CWindowingTools(); + const su2double weight = windowEvaluator.GetWndWeight(config->GetKindWindow(), TimeIter, IterAvg_Obj - 1); + seeding = weight / IterAvg_Obj; } else { seeding = 0.0; } } - - if (rank == MASTER_NODE){ + if (rank == MASTER_NODE) { SU2_TYPE::SetDerivative(ObjFunc, SU2_TYPE::GetValue(seeding)); } else { SU2_TYPE::SetDerivative(ObjFunc, 0.0); diff --git a/SU2_CFD/src/integration/CStructuralIntegration.cpp b/SU2_CFD/src/integration/CStructuralIntegration.cpp index 1dd5573923cd..5fc18babb3f2 100644 --- a/SU2_CFD/src/integration/CStructuralIntegration.cpp +++ b/SU2_CFD/src/integration/CStructuralIntegration.cpp @@ -134,9 +134,6 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * /*--- Set the Jacobian according to the different time integration methods ---*/ switch (config->GetKind_TimeIntScheme_FEA()) { - case (STRUCT_TIME_INT::CD_EXPLICIT): - solver_container[MainSolver]->ImplicitNewmark_Iteration(geometry, numerics, config); - break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): solver_container[MainSolver]->ImplicitNewmark_Iteration(geometry, numerics, config); break; @@ -161,16 +158,13 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * } } - /*--- Solver linear system ---*/ + /*--- Solve linear system ---*/ solver_container[MainSolver]->Solve_System(geometry, config); /*--- Update solution ---*/ switch (config->GetKind_TimeIntScheme_FEA()) { - case (STRUCT_TIME_INT::CD_EXPLICIT): - solver_container[MainSolver]->ImplicitNewmark_Update(geometry, config); - break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): solver_container[MainSolver]->ImplicitNewmark_Update(geometry, config); break; @@ -193,15 +187,13 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * void CStructuralIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolver *solver, const CConfig *config, unsigned short iMesh) { - bool fsi = config->GetFSI_Simulation(); + const bool fsi = config->GetFSI_Simulation(); /*--- Update the solution according to the integration scheme used ---*/ switch (config->GetKind_TimeIntScheme_FEA()) { - case (STRUCT_TIME_INT::CD_EXPLICIT): - break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): - if (fsi) solver->ImplicitNewmark_Relaxation(geometry, config); + if (fsi && config->GetRelaxation()) solver->ImplicitNewmark_Relaxation(geometry, config); break; case (STRUCT_TIME_INT::GENERALIZED_ALPHA): solver->GeneralizedAlpha_UpdateSolution(geometry, config); @@ -215,12 +207,5 @@ void CStructuralIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolv /*--- If FSI problem, save the last Aitken relaxation parameter of the previous time step ---*/ - if (fsi) { - - su2double WAitk=0.0; - - WAitk = solver->GetWAitken_Dyn(); - solver->SetWAitken_Dyn_tn1(WAitk); - - } + if (fsi) solver->SetWAitken_Dyn_tn1(); } diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index d95bb5a98415..1aaa0d422fe5 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -132,7 +132,7 @@ void CFlowTractionInterface::GetPhysical_Constants(CSolver *flow_solution, CSolv /*--- Apply a ramp to the transfer of the fluid loads ---*/ su2double ModAmpl = 0.0; - su2double CurrentTime = struct_config->GetCurrent_DynTime(); + su2double CurrentTime = struct_config->GetCurrent_UnstTime(); bool Ramp_Load = struct_config->GetRamp_Load(); su2double Ramp_Time = struct_config->GetRamp_Time(); diff --git a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp index 4a6e61dfc2c9..a37971f0aeb8 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp @@ -32,34 +32,6 @@ CDiscAdjFEAIteration::CDiscAdjFEAIteration(const CConfig *config) : CIteration(config), CurrentRecording(NONE) { fem_iteration = new CFEAIteration(config); - - // TEMPORARY output only for standalone structural problems - if ((!config->GetFSI_Simulation()) && (rank == MASTER_NODE)) { - bool de_effects = config->GetDE_Effects(); - unsigned short iVar; - - /*--- Header of the temporary output file ---*/ - ofstream myfile_res; - myfile_res.open("Results_Reverse_Adjoint.txt"); - - myfile_res << "Obj_Func" - << " "; - for (iVar = 0; iVar < config->GetnElasticityMod(); iVar++) myfile_res << "Sens_E_" << iVar << "\t"; - - for (iVar = 0; iVar < config->GetnPoissonRatio(); iVar++) myfile_res << "Sens_Nu_" << iVar << "\t"; - - if (config->GetTime_Domain()) { - for (iVar = 0; iVar < config->GetnMaterialDensity(); iVar++) myfile_res << "Sens_Rho_" << iVar << "\t"; - } - - if (de_effects) { - for (iVar = 0; iVar < config->GetnElectric_Field(); iVar++) myfile_res << "Sens_EField_" << iVar << "\t"; - } - - myfile_res << endl; - - myfile_res.close(); - } } CDiscAdjFEAIteration::~CDiscAdjFEAIteration() = default; @@ -68,44 +40,54 @@ void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integrat CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short iZone, unsigned short iInst) { - unsigned long iPoint; auto solvers0 = solver[iZone][iInst][MESH_0]; auto geometry0 = geometry[iZone][iInst][MESH_0]; auto dirNodes = solvers0[FEA_SOL]->GetNodes(); auto adjNodes = solvers0[ADJFEA_SOL]->GetNodes(); - /*--- For the dynamic adjoint, load direct solutions from restart files. ---*/ + auto StoreDirectSolution = [&]() { + for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) { + adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution(iPoint)); + } + }; + + /*--- For the dynamic adjoint, load direct solutions from restart files. + * For steady, store the direct solution to be able to reset it later. ---*/ if (config[iZone]->GetTime_Domain()) { const int TimeIter = config[iZone]->GetTimeIter(); const int Direct_Iter = SU2_TYPE::Int(config[iZone]->GetUnst_AdjointIter()) - TimeIter - 1; - /*--- We want to load the already converged solution at timesteps n and n-1 ---*/ - - /*--- Load solution at timestep n-1 ---*/ + /*--- We need the already converged solution at timesteps n and n-1. On the first + * adjoint time step we load both, then n-1 becomes "n" and we only load n-2. ---*/ - LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 1); + if (TimeIter > 0) { + /*--- Save n-1 to become n. ---*/ + for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) { + adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution_time_n(iPoint)); + } + } - /*--- Push solution back to correct array ---*/ + /*--- Load solution at timestep n-1 and push back to correct array. ---*/ + LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter - 1); dirNodes->Set_Solution_time_n(); - /*--- Load solution timestep n ---*/ + /*--- Load or set solution at timestep n. ---*/ - LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter); - - /*--- Store FEA solution also in the adjoint solver in order to be able to reset it later ---*/ - - for (iPoint = 0; iPoint < geometry0->GetnPoint(); iPoint++) { - adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution(iPoint)); + if (TimeIter == 0) { + LoadDynamic_Solution(geometry, solver, config, iZone, iInst, Direct_Iter); + StoreDirectSolution(); + } else { + /*--- Set n-1 as n. ---*/ + for (auto iPoint = 0ul; iPoint < geometry0->GetnPoint(); iPoint++) + for (auto iVar = 0u; iVar < solvers0[ADJFEA_SOL]->GetnVar(); iVar++) + dirNodes->SetSolution(iPoint, iVar, adjNodes->GetSolution_Direct(iPoint)[iVar]); } } else { - /*--- Store FEA solution also in the adjoint solver in order to be able to reset it later ---*/ - - for (iPoint = 0; iPoint < geometry0->GetnPoint(); iPoint++) { - adjNodes->SetSolution_Direct(iPoint, dirNodes->GetSolution(iPoint)); - } + /*--- Steady. ---*/ + StoreDirectSolution(); } solvers0[ADJFEA_SOL]->Preprocessing(geometry0, solvers0, config[iZone], MESH_0, 0, RUNTIME_ADJFEA_SYS, false); @@ -115,27 +97,19 @@ void CDiscAdjFEAIteration::Preprocess(COutput* output, CIntegration**** integrat void CDiscAdjFEAIteration::LoadDynamic_Solution(CGeometry**** geometry, CSolver***** solver, CConfig** config, unsigned short iZone, unsigned short iInst, int val_DirectIter) { - unsigned short iVar; - unsigned long iPoint; - bool update_geo = false; // TODO: check + /*--- Set to false to prevent updating Solution_time_n when loading primal solutions of unsteady cases. ---*/ + const bool update_geo = false; + auto*** solvers = solver[iZone][iInst]; if (val_DirectIter >= 0) { if (rank == MASTER_NODE && iZone == ZONE_0) - cout << " Loading FEA solution from direct iteration " << val_DirectIter << "." << endl; - solver[iZone][iInst][MESH_0][FEA_SOL]->LoadRestart( - geometry[iZone][iInst], solver[iZone][iInst], config[iZone], val_DirectIter, update_geo); + cout << " Loading FEA solution from direct iteration " << val_DirectIter << ".\n"; + solvers[MESH_0][FEA_SOL]->LoadRestart(geometry[iZone][iInst], solvers, config[iZone], val_DirectIter, update_geo); } else { - /*--- If there is no solution file we set the freestream condition ---*/ + /*--- If there is no solution file we set the initial conditions. ---*/ if (rank == MASTER_NODE && iZone == ZONE_0) - cout << " Setting static conditions at direct iteration " << val_DirectIter << "." << endl; - /*--- Push solution back to correct array ---*/ - for (iPoint = 0; iPoint < geometry[iZone][iInst][MESH_0]->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < solver[iZone][iInst][MESH_0][FEA_SOL]->GetnVar(); iVar++) { - solver[iZone][iInst][MESH_0][FEA_SOL]->GetNodes()->SetSolution(iPoint, iVar, 0.0); - solver[iZone][iInst][MESH_0][FEA_SOL]->GetNodes()->SetSolution_Accel(iPoint, iVar, 0.0); - solver[iZone][iInst][MESH_0][FEA_SOL]->GetNodes()->SetSolution_Vel(iPoint, iVar, 0.0); - } - } + cout << " Setting static conditions at direct iteration " << val_DirectIter << ".\n"; + solvers[MESH_0][FEA_SOL]->SetInitialCondition(geometry[iZone][iInst], solvers, config[iZone], val_DirectIter); } } @@ -196,7 +170,7 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** ge const int mat_knowles = MAT_KNOWLES+offset; const int de_term = DE_TERM+offset; - for (unsigned short iProp = 0; iProp < config[iZone]->GetnElasticityMod(); iProp++) { + for (unsigned short iProp = 0; iProp < config[iZone]->GetnElasticityMat(); iProp++) { su2double E = adj_solver->GetVal_Young(iProp); su2double nu = adj_solver->GetVal_Poisson(iProp); su2double rho = adj_solver->GetVal_Rho(iProp); @@ -320,70 +294,27 @@ void CDiscAdjFEAIteration::Postprocess(COutput* output, CIntegration**** integra CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short iZone, unsigned short iInst) { - const bool dynamic = (config[iZone]->GetTime_Domain()); auto solvers0 = solver[iZone][iInst][MESH_0]; - // TEMPORARY output only for standalone structural problems - if (config[iZone]->GetAdvanced_FEAElementBased() && (rank == MASTER_NODE)) { - unsigned short iVar; - - const bool de_effects = config[iZone]->GetDE_Effects(); - - /*--- Header of the temporary output file ---*/ - ofstream myfile_res; - myfile_res.open("Results_Reverse_Adjoint.txt", ios::app); - - myfile_res.precision(15); - - myfile_res << config[iZone]->GetTimeIter() << "\t"; - - solvers0[FEA_SOL]->Evaluate_ObjFunc(config[iZone], solvers0); - myfile_res << scientific << solvers0[FEA_SOL]->GetTotal_ComboObj() << "\t"; - - for (iVar = 0; iVar < config[iZone]->GetnElasticityMod(); iVar++) - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_E(iVar) << "\t"; - for (iVar = 0; iVar < config[iZone]->GetnPoissonRatio(); iVar++) - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_Nu(iVar) << "\t"; - if (dynamic) { - for (iVar = 0; iVar < config[iZone]->GetnMaterialDensity(); iVar++) - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_Rho(iVar) << "\t"; - } - if (de_effects) { - for (iVar = 0; iVar < config[iZone]->GetnElectric_Field(); iVar++) - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_EField(iVar) << "\t"; - } - for (iVar = 0; iVar < solvers0[ADJFEA_SOL]->GetnDVFEA(); iVar++) { - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_DVFEA(iVar) << "\t"; - } - - myfile_res << endl; - - myfile_res.close(); - } - // TEST: for implementation of python framework in standalone structural problems if (config[iZone]->GetAdvanced_FEAElementBased() && (rank == MASTER_NODE)) { /*--- Header of the temporary output file ---*/ ofstream myfile_res; - bool outputDVFEA = false; + bool outputDVFEA = true; switch (config[iZone]->GetDV_FEA()) { case YOUNG_MODULUS: myfile_res.open("grad_young.opt"); - outputDVFEA = true; break; case POISSON_RATIO: myfile_res.open("grad_poisson.opt"); - outputDVFEA = true; break; case DENSITY_VAL: case DEAD_WEIGHT: myfile_res.open("grad_density.opt"); - outputDVFEA = true; break; case ELECTRIC_FIELD: myfile_res.open("grad_efield.opt"); - outputDVFEA = true; break; default: outputDVFEA = false; @@ -391,23 +322,14 @@ void CDiscAdjFEAIteration::Postprocess(COutput* output, CIntegration**** integra } if (outputDVFEA) { - unsigned short iDV; - unsigned short nDV = solvers0[ADJFEA_SOL]->GetnDVFEA(); - - myfile_res << "INDEX" - << "\t" - << "GRAD" << endl; - + const auto nDV = solvers0[ADJFEA_SOL]->GetnDVFEA(); + myfile_res << "INDEX\tGRAD\n"; myfile_res.precision(15); - for (iDV = 0; iDV < nDV; iDV++) { - myfile_res << iDV; - myfile_res << "\t"; - myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_DVFEA(iDV); - myfile_res << endl; + for (auto iDV = 0u; iDV < nDV; iDV++) { + myfile_res << iDV << "\t"; + myfile_res << scientific << solvers0[ADJFEA_SOL]->GetTotal_Sens_DVFEA(iDV) << "\n"; } - - myfile_res.close(); } } } diff --git a/SU2_CFD/src/iteration/CFEAIteration.cpp b/SU2_CFD/src/iteration/CFEAIteration.cpp index 0190d2d52347..91b9e7156813 100644 --- a/SU2_CFD/src/iteration/CFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CFEAIteration.cpp @@ -69,8 +69,8 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom } else if (nonlinear && !incremental_load) { /*--- THIS IS THE DIRECT APPROACH (NO INCREMENTAL LOAD APPLIED) ---*/ - /*--- Keep the current inner iter, we need to restore it in discrete adjoint cases as file output depends on it - * ---*/ + /*--- Keep the current inner iter, we need to restore it in discrete adjoint cases + * because file output depends on it. ---*/ const auto CurIter = config[val_iZone]->GetInnerIter(); /*--- Newton-Raphson subiterations ---*/ @@ -84,11 +84,11 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom if (disc_adj_fem) { config[val_iZone]->SetInnerIter(CurIter); break; - } StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, - FFDBox, val_iZone, INST_0); + } + StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, + FFDBox, val_iZone, INST_0); - if (StopCalc && (IntIter > 0)) break; - + if (StopCalc && (IntIter > 0)) break; } } else { @@ -101,13 +101,16 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom /*--- Run two nonlinear iterations to check if incremental loading can be skipped ---*/ - for (IntIter = 0; IntIter < 2; ++IntIter) { + auto Iterate = [&](unsigned long IntIter) { config[val_iZone]->SetInnerIter(IntIter); - feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, val_iInst); StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox, val_iZone, INST_0); + }; + + for (IntIter = 0; IntIter < 2; ++IntIter) { + Iterate(IntIter); } /*--- Early return if we already meet the convergence criteria. ---*/ @@ -125,13 +128,7 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom /*--- Newton-Raphson subiterations ---*/ for (IntIter = 2; IntIter < config[val_iZone]->GetnInner_Iter(); IntIter++) { - config[val_iZone]->SetInnerIter(IntIter); - - feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, val_iInst); - - StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, - FFDBox, val_iZone, INST_0); - + Iterate(IntIter); if (StopCalc) break; } @@ -162,14 +159,7 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeom /*--- Newton-Raphson subiterations ---*/ for (IntIter = 0; IntIter < config[val_iZone]->GetnInner_Iter(); IntIter++) { - config[val_iZone]->SetInnerIter(IntIter); - - feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, - val_iInst); - - StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, - FFDBox, val_iZone, INST_0); - + Iterate(IntIter); if (StopCalc && (IntIter > 0)) break; } } @@ -183,21 +173,16 @@ void CFEAIteration::Update(COutput* output, CIntegration**** integration, CGeome CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) { - const bool dynamic = (config[val_iZone]->GetTime_Domain()); - const bool fsi = config[val_iZone]->GetFSI_Simulation(); - CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL]; /*----------------- Update structural solver ----------------------*/ - if (dynamic) { - integration[val_iZone][val_iInst][FEA_SOL]->SetDualTime_Solver( - geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0][FEA_SOL], config[val_iZone], - MESH_0); - - } else if (fsi) { - /*--- For FSI problems, output the relaxed result, which is the one transferred into the fluid domain (for restart - * purposes) ---*/ + if (config[val_iZone]->GetTime_Domain()) { + integration[val_iZone][val_iInst][FEA_SOL]->SetDualTime_Solver(geometry[val_iZone][val_iInst][MESH_0], feaSolver, + config[val_iZone], MESH_0); + } else if (config[val_iZone]->GetFSI_Simulation() && config[val_iZone]->GetRelaxation()) { + /*--- For FSI problems with relaxation, output the relaxed result, which is the one transferred into the fluid + * domain (for consistent restart purposes). ---*/ if (config[val_iZone]->GetKind_TimeIntScheme_FEA() == STRUCT_TIME_INT::NEWMARK_IMPLICIT) { feaSolver->ImplicitNewmark_Relaxation(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); } diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp index 6fd3a35042c1..fc87cbc19a58 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp @@ -38,23 +38,19 @@ CFEAElasticity::CFEAElasticity(unsigned short val_nDim, unsigned short val_nVar, bool body_forces = config->GetDeadLoad(); // Body forces (dead loads). bool pseudo_static = config->GetPseudoStatic(); - unsigned short iVar, nProp; + unsigned short iVar; /*--- Initialize vector structures for multiple material definition ---*/ - nProp = config->GetnElasticityMod(); + const auto nProp = config->GetnElasticityMat(); E_i = new su2double[nProp]; for (iVar = 0; iVar < nProp; iVar++) E_i[iVar] = config->GetElasticyMod(iVar); - nProp = config->GetnPoissonRatio(); - Nu_i = new su2double[nProp]; for (iVar = 0; iVar < nProp; iVar++) Nu_i[iVar] = config->GetPoissonRatio(iVar); - nProp = config->GetnMaterialDensity(); - Rho_s_i = new su2double[nProp]; // For inertial effects Rho_s_DL_i = new su2double[nProp]; // For dead loads diff --git a/SU2_CFD/src/output/CAdjElasticityOutput.cpp b/SU2_CFD/src/output/CAdjElasticityOutput.cpp index 8ed7017aba2e..e30f2895f654 100644 --- a/SU2_CFD/src/output/CAdjElasticityOutput.cpp +++ b/SU2_CFD/src/output/CAdjElasticityOutput.cpp @@ -27,6 +27,7 @@ #include "../../include/output/CAdjElasticityOutput.hpp" +#include #include "../../../Common/include/geometry/CGeometry.hpp" #include "../../include/solvers/CSolver.hpp" @@ -50,8 +51,8 @@ CAdjElasticityOutput::CAdjElasticityOutput(CConfig *config, unsigned short nDim) requestedScreenFields.emplace_back("INNER_ITER"); requestedScreenFields.emplace_back("ADJOINT_DISP_X"); requestedScreenFields.emplace_back("ADJOINT_DISP_Y"); - requestedScreenFields.emplace_back("SENS_E"); - requestedScreenFields.emplace_back("SENS_NU"); + requestedScreenFields.emplace_back("SENS_E_0"); + requestedScreenFields.emplace_back("SENS_NU_0"); nRequestedScreenFields = requestedScreenFields.size(); } @@ -62,9 +63,9 @@ CAdjElasticityOutput::CAdjElasticityOutput(CConfig *config, unsigned short nDim) nRequestedVolumeFields = requestedVolumeFields.size(); } - if (find(requestedVolumeFields.begin(), requestedVolumeFields.end(), string("SENSITIVITY")) == requestedVolumeFields.end()) { + if (find(requestedVolumeFields.begin(), requestedVolumeFields.end(), "SENSITIVITY") == requestedVolumeFields.end()) { requestedVolumeFields.emplace_back("SENSITIVITY"); - nRequestedVolumeFields ++; + nRequestedVolumeFields++; } stringstream ss; @@ -97,47 +98,69 @@ CAdjElasticityOutput::~CAdjElasticityOutput() = default; void CAdjElasticityOutput::SetHistoryOutputFields(CConfig *config){ - // Residuals + /*--- Residuals ---*/ AddHistoryOutput("ADJOINT_DISP_X", "rms[Ux_adj]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint of the X displacements.", HistoryFieldType::RESIDUAL); AddHistoryOutput("ADJOINT_DISP_Y", "rms[Uy_adj]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint of the Y displacements.", HistoryFieldType::RESIDUAL); - AddHistoryOutput("ADJOINT_DISP_Z", "rms[Uz_adj]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint of the Z displacements.", HistoryFieldType::RESIDUAL); + if (nVar_FEM == 3) { + AddHistoryOutput("ADJOINT_DISP_Z", "rms[Uz_adj]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint of the Z displacements.", HistoryFieldType::RESIDUAL); + } - //Sensitivities - AddHistoryOutput("SENS_E", "Sens[E]", ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Elasticity modulus"); - AddHistoryOutput("SENS_NU", "Sens[Nu]", ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Poisson ratio"); + /*--- Sensitivities ---*/ + for (auto iVar = 0u; iVar < config->GetnElasticityMat(); iVar++) { + const auto iVarS = std::to_string(iVar); + AddHistoryOutput("SENS_E_" + iVarS, "Sens[E" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Elasticity modulus"); + AddHistoryOutput("SENS_NU_" + iVarS, "Sens[Nu" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Poisson ratio"); + if (config->GetTime_Domain() && !config->GetPseudoStatic()) { + AddHistoryOutput("SENS_RHO_" + iVarS, "Sens[Rho" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Material density"); + } + if (config->GetDeadLoad()) { + AddHistoryOutput("SENS_RHO_DL_" + iVarS, "Sens[RhoDL" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Dead load density"); + } + } + if (config->GetDE_Effects()) { + for (auto iVar = 0u; iVar < config->GetnElectric_Field(); iVar++) { + const auto iVarS = std::to_string(iVar); + AddHistoryOutput("SENS_EFIELD_" + iVarS, "Sens[EField" + iVarS + ']', ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", "d Objective / d Electric field"); + } + } AddHistoryOutput("LINSOL_ITER", "LinSolIter", ScreenOutputFormat::INTEGER, "LINSOL", "Number of iterations of the linear solver."); AddHistoryOutput("LINSOL_RESIDUAL", "LinSolRes", ScreenOutputFormat::FIXED, "LINSOL", "Residual of the linear solver."); - AddHistoryOutput("BGS_ADJ_DISP_X", "bgs[A_Ux]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint X displacement.", HistoryFieldType::RESIDUAL); - AddHistoryOutput("BGS_ADJ_DISP_Y", "bgs[A_Uy]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Y displacement.", HistoryFieldType::RESIDUAL); - AddHistoryOutput("BGS_ADJ_DISP_Z", "bgs[A_Uz]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Z displacement.", HistoryFieldType::RESIDUAL); - + if (multiZone) { + AddHistoryOutput("BGS_ADJ_DISP_X", "bgs[A_Ux]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint X displacement.", HistoryFieldType::RESIDUAL); + AddHistoryOutput("BGS_ADJ_DISP_Y", "bgs[A_Uy]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Y displacement.", HistoryFieldType::RESIDUAL); + if (nVar_FEM == 3) { + AddHistoryOutput("BGS_ADJ_DISP_Z", "bgs[A_Uz]", ScreenOutputFormat::FIXED, "BGS_RES", "BGS residual of the adjoint Z displacement.", HistoryFieldType::RESIDUAL); + } + } } inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) { SetHistoryOutputValue("ADJOINT_DISP_X", log10(solver[ADJFEA_SOL]->GetRes_RMS(0))); SetHistoryOutputValue("ADJOINT_DISP_Y", log10(solver[ADJFEA_SOL]->GetRes_RMS(1))); - if (nVar_FEM == 3){ + if (nVar_FEM == 3) { SetHistoryOutputValue("ADJOINT_DISP_Z", log10(solver[ADJFEA_SOL]->GetRes_RMS(2))); } - su2double Total_SensE = 0.0; su2double Total_SensNu = 0.0; - if (config->GetnElasticityMod() == 1) { - Total_SensE = solver[ADJFEA_SOL]->GetGlobal_Sens_E(0); - Total_SensNu = solver[ADJFEA_SOL]->GetGlobal_Sens_Nu(0); + for (unsigned short iVar = 0; iVar < config->GetnElasticityMat(); iVar++) { + const auto iVarS = std::to_string(iVar); + SetHistoryOutputValue("SENS_E_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_E(iVar)); + SetHistoryOutputValue("SENS_NU_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_Nu(iVar)); + if (config->GetTime_Domain() && !config->GetPseudoStatic()) { + SetHistoryOutputValue("SENS_RHO_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_Rho(iVar)); + } + if (config->GetDeadLoad()) { + SetHistoryOutputValue("SENS_RHO_DL_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_Rho_DL(iVar)); + } } - else { - for (unsigned short iVar = 0; iVar < config->GetnElasticityMod(); iVar++){ - Total_SensE += pow(solver[ADJFEA_SOL]->GetGlobal_Sens_E(iVar),2); - Total_SensNu += pow(solver[ADJFEA_SOL]->GetGlobal_Sens_Nu(iVar),2); + if (config->GetDE_Effects()) { + for (auto iVar = 0u; iVar < config->GetnElectric_Field(); iVar++) { + const auto iVarS = std::to_string(iVar); + SetHistoryOutputValue("SENS_EFIELD_" + iVarS, solver[ADJFEA_SOL]->GetTotal_Sens_EField(iVar)); } - Total_SensE = sqrt(Total_SensE); - Total_SensNu = sqrt(Total_SensNu); } - SetHistoryOutputValue("SENS_E", Total_SensE); - SetHistoryOutputValue("SENS_NU", Total_SensNu); SetHistoryOutputValue("LINSOL_ITER", solver[ADJFEA_SOL]->GetIterLinSolver()); SetHistoryOutputValue("LINSOL_RESIDUAL", log10(solver[ADJFEA_SOL]->GetResLinSolver())); @@ -170,6 +193,23 @@ void CAdjElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, SetVolumeOutputValue("SENSITIVITY-Y", iPoint, Node_Struc->GetSensitivity(iPoint, 1)); if (nDim == 3) SetVolumeOutputValue("SENSITIVITY-Z", iPoint, Node_Struc->GetSensitivity(iPoint, 2)); + + if (!config->GetTime_Domain()) return; + + SetVolumeOutputValue("SENS_DISP-X", iPoint, Node_Struc->GetSolution_time_n(iPoint, 0)); + SetVolumeOutputValue("SENS_DISP-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, 1)); + if (nDim == 3) + SetVolumeOutputValue("SENS_DISP-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2)); + + SetVolumeOutputValue("SENS_VEL-X", iPoint, Node_Struc->GetSolution_time_n(iPoint, nDim)); + SetVolumeOutputValue("SENS_VEL-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, nDim + 1)); + if (nDim == 3) + SetVolumeOutputValue("SENS_VEL-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 5)); + + SetVolumeOutputValue("SENS_ACCEL-X", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2 * nDim)); + SetVolumeOutputValue("SENS_ACCEL-Y", iPoint, Node_Struc->GetSolution_time_n(iPoint, 2 * nDim + 1)); + if (nDim == 3) + SetVolumeOutputValue("SENS_ACCEL-Z", iPoint, Node_Struc->GetSolution_time_n(iPoint, 8)); } void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){ @@ -200,4 +240,23 @@ void CAdjElasticityOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("SENSITIVITY-Z", "Sensitivity_z", "SENSITIVITY", "geometric sensitivity in the z direction"); /// END_GROUP + if (!config->GetTime_Domain()) return; + + /*--- Sensitivities with respect to initial conditions. ---*/ + + AddVolumeOutput("SENS_DISP-X", "SensInitialDisp_x", "SENSITIVITY_T0", "sensitivity to the initial x displacement"); + AddVolumeOutput("SENS_DISP-Y", "SensInitialDisp_y", "SENSITIVITY_T0", "sensitivity to the initial y displacement"); + if (nDim == 3) + AddVolumeOutput("SENS_DISP-Z", "SensInitialDisp_z", "SENSITIVITY_T0", "sensitivity to the initial z displacement"); + + AddVolumeOutput("SENS_VEL-X", "SensInitialVel_x", "SENSITIVITY_T0", "sensitivity to the initial x velocity"); + AddVolumeOutput("SENS_VEL-Y", "SensInitialVel_y", "SENSITIVITY_T0", "sensitivity to the initial y velocity"); + if (nDim == 3) + AddVolumeOutput("SENS_VEL-Z", "SensInitialVel_z", "SENSITIVITY_T0", "sensitivity to the initial z velocity"); + + AddVolumeOutput("SENS_ACCEL-X", "SensInitialAccel_x", "SENSITIVITY_T0", "sensitivity to the initial x acceleration"); + AddVolumeOutput("SENS_ACCEL-Y", "SensInitialAccel_y", "SENSITIVITY_T0", "sensitivity to the initial y acceleration"); + if (nDim == 3) + AddVolumeOutput("SENS_ACCEL-Z", "SensInitialAccel_z", "SENSITIVITY_T0", "sensitivity to the initial z acceleration"); + } diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index 3fc02d126e75..bff61c408c9f 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -873,8 +873,8 @@ void CFlowOutput::SetCustomOutputs(const CSolver* const* solver, const CGeometry const auto varIdx = i % CustomOutput::MAX_VARS_PER_SOLVER; if (solIdx == FLOW_SOL) { return flowNodes->GetPrimitive(iPoint, varIdx); - } return solver[solIdx]->GetNodes()->GetSolution(iPoint, varIdx); - + } + return solver[solIdx]->GetNodes()->GetSolution(iPoint, varIdx); } else { return *output.otherOutputs[i - CustomOutput::NOT_A_VARIABLE]; } diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index f11c9f580b66..b64d409e5cc5 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -84,16 +84,7 @@ CDiscAdjFEASolver::CDiscAdjFEASolver(CGeometry *geometry, CConfig *config, CSolv /*--- Initialize vector structures for multiple material definition ---*/ - nMPROP = config->GetnElasticityMod(); - - /*--- For a material to be fully defined, we need to have the same number for all three parameters ---*/ - bool checkDef = ((config->GetnElasticityMod() == config->GetnPoissonRatio()) && - (config->GetnElasticityMod() == config->GetnMaterialDensity()) && - (config->GetnMaterialDensity() == config->GetnPoissonRatio())); - - if (!checkDef){ - SU2_MPI::Error("WARNING: For a material to be fully defined, E, Nu and Rho need to have the same dimensions.", CURRENT_FUNCTION); - } + nMPROP = config->GetnElasticityMat(); E.resize(nMPROP); Nu.resize(nMPROP); @@ -133,28 +124,21 @@ CDiscAdjFEASolver::~CDiscAdjFEASolver() { delete nodes; } void CDiscAdjFEASolver::SetRecording(CGeometry* geometry, CConfig *config){ - unsigned long iPoint; - unsigned short iVar; - /*--- Reset the solution to the initial (converged) solution ---*/ - for (iPoint = 0; iPoint < nPoint; iPoint++) - for (iVar = 0; iVar < nVar; iVar++) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + for (auto iVar = 0u; iVar < nVar; iVar++) direct_solver->GetNodes()->SetSolution(iPoint, iVar, nodes->GetSolution_Direct(iPoint)[iVar]); + } /*--- Reset the input for time n ---*/ if (config->GetTime_Domain()) { - for (iPoint = 0; iPoint < nPoint; iPoint++) - for (iVar = 0; iVar < nVar; iVar++) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) + for (auto iVar = 0u; iVar < nVar; iVar++) AD::ResetInput(direct_solver->GetNodes()->GetSolution_time_n(iPoint)[iVar]); } - /*--- Set the Jacobian to zero since this is not done inside the meanflow iteration - * when running the discrete adjoint solver. ---*/ - - direct_solver->Jacobian.SetValZero(); - /*--- Set indices to zero ---*/ RegisterVariables(geometry, config, true); @@ -214,10 +198,11 @@ void CDiscAdjFEASolver::RegisterVariables(CGeometry *geometry, CConfig *config, Rho_DL.Register(); if (de_effects) EField.Register(); if (fea_dv) DV.Register(); + } - /*--- Register the flow tractions ---*/ - if (config->GetnMarker_Fluid_Load() > 0) - direct_solver->GetNodes()->RegisterFlowTraction(); + /*--- Register or reset the flow tractions ---*/ + if (config->GetnMarker_Fluid_Load() > 0) { + direct_solver->GetNodes()->RegisterFlowTraction(reset); } } @@ -303,7 +288,7 @@ void CDiscAdjFEASolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *c if (config->GetnMarker_Fluid_Load() > 0) { for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){ for (unsigned short iDim = 0; iDim < nDim; iDim++){ - su2double val_sens = direct_solver->GetNodes()->ExtractFlowTraction_Sensitivity(iPoint,iDim); + su2double val_sens = direct_solver->GetNodes()->ExtractFlowTractionSensitivity(iPoint,iDim); nodes->SetFlowTractionSensitivity(iPoint, iDim, val_sens); } } @@ -355,12 +340,17 @@ void CDiscAdjFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_cont void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*){ - E.UpdateTotal(); - Nu.UpdateTotal(); - Rho.UpdateTotal(); - Rho_DL.UpdateTotal(); - if (de_effects) EField.UpdateTotal(); - if (fea_dv) DV.UpdateTotal(); + const bool time_domain = config->GetTime_Domain(); + + /*--- Store the final material sensitivities for the time step to increment them on the next time step. ---*/ + if (time_domain) { + E.Store(); + Nu.Store(); + Rho.Store(); + Rho_DL.Store(); + if (de_effects) EField.Store(); + if (fea_dv) DV.Store(); + } /*--- Extract the topology optimization density sensitivities. ---*/ @@ -368,8 +358,6 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo /*--- Extract the geometric sensitivities ---*/ - const bool time_domain = config->GetTime_Domain(); - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { auto Coord = geometry->nodes->GetCoord(iPoint); diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index d6d94cf0d58b..e236b19114d1 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -146,11 +146,6 @@ void CDiscAdjSolver::SetRecording(CGeometry* geometry, CConfig *config){ END_SU2_OMP_FOR } - /*--- Set the Jacobian to zero since this is not done inside the fluid iteration - * when running the discrete adjoint solver. ---*/ - - direct_solver->Jacobian.SetValZero(); - /*--- Set indices to zero ---*/ RegisterVariables(geometry, config, true); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 57495c3d8297..858257f094c9 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -50,6 +50,7 @@ CFEASolver::CFEASolver(LINEAR_SOLVER_MODE mesh_deform_mode) : CFEASolverBase(mes CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CFEASolverBase(geometry, config) { bool dynamic = (config->GetTime_Domain()); + config->SetDelta_UnstTimeND(config->GetDelta_UnstTime()); /*--- Test whether we consider dielectric elastomers ---*/ bool de_effects = config->GetDE_Effects(); @@ -382,8 +383,8 @@ void CFEASolver::Set_ElementProperties(CGeometry *geometry, CConfig *config) { /*--- Detect a wrong solution file ---*/ if (iElem_Global_Local != nElement) { - SU2_MPI::Error(string("The properties file ") + filename + string(" doesn't match with the mesh file!\n") + - string("It could be empty lines at the end of the file."), CURRENT_FUNCTION); + SU2_MPI::Error("The properties file " + filename + " doesn't match with the mesh file!\n" + "It could be empty lines at the end of the file.", CURRENT_FUNCTION); } } @@ -614,8 +615,8 @@ void CFEASolver::SetInitialCondition(CGeometry **geometry, CSolver ***solver_con SU2_OMP_PARALLEL { + su2double zeros[MAXNVAR] = {0.0}; if (!config->GetPrestretch()) { - su2double zeros[MAXNVAR] = {0.0}; SU2_OMP_FOR_STAT(omp_chunk_size) for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) nodes->SetSolution(iPoint, zeros); @@ -627,6 +628,14 @@ void CFEASolver::SetInitialCondition(CGeometry **geometry, CSolver ***solver_con nodes->SetSolution(iPoint, nodes->GetPrestretch(iPoint)); END_SU2_OMP_FOR } + if (config->GetTime_Domain()) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { + nodes->SetSolution_Vel(iPoint, zeros); + nodes->SetSolution_Accel(iPoint, zeros); + } + END_SU2_OMP_FOR + } } END_SU2_OMP_PARALLEL } @@ -1282,7 +1291,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, if (outputReactions) { - bool dynamic = (config->GetDynamic_Analysis() == DYNAMIC); + const bool dynamic = config->GetTime_Domain(); ofstream myfile; myfile.open ("Reactions.txt"); @@ -1331,9 +1340,6 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, else if (dynamic) { switch (config->GetKind_TimeIntScheme_FEA()) { - case (STRUCT_TIME_INT::CD_EXPLICIT): - cout << "NOT IMPLEMENTED YET" << endl; - break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): /*--- Loop over all points, and set aux vector TimeRes_Aux = a0*U+a2*U'+a3*U'' ---*/ @@ -1481,14 +1487,11 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, con void CFEASolver::Compute_IntegrationConstants(const CConfig *config) { - su2double Delta_t= config->GetDelta_DynTime(); + su2double Delta_t= config->GetDelta_UnstTime(); su2double gamma = config->GetNewmark_gamma(), beta = config->GetNewmark_beta(); switch (config->GetKind_TimeIntScheme_FEA()) { - case (STRUCT_TIME_INT::CD_EXPLICIT): - cout << "NOT IMPLEMENTED YET" << endl; - break; case (STRUCT_TIME_INT::NEWMARK_IMPLICIT): /*--- Integration constants for Newmark scheme ---*/ @@ -1652,7 +1655,7 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, const CConfig *config, unsigned const su2double *DispDirLocal = config->GetDisp_Dir(TagBound); su2double DispDirMod = Norm(nDim, DispDirLocal); - su2double CurrentTime = config->GetCurrent_DynTime(); + su2double CurrentTime = config->GetCurrent_UnstTime(); su2double RampTime = config->GetRamp_Time(); su2double ModAmpl = Compute_LoadCoefficient(CurrentTime, RampTime, config); @@ -1818,7 +1821,7 @@ void CFEASolver::BC_Normal_Load(CGeometry *geometry, const CConfig *config, unsi /*--- Retrieve the normal pressure and the application conditions for the considered boundary. ---*/ - su2double CurrentTime = config->GetCurrent_DynTime(); + su2double CurrentTime = config->GetCurrent_UnstTime(); su2double Ramp_Time = config->GetRamp_Time(); su2double ModAmpl = Compute_LoadCoefficient(CurrentTime, Ramp_Time, config); @@ -1914,7 +1917,7 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, const CConfig *config, unsigne /*--- Compute the norm of the vector that was passed in the config file. ---*/ su2double LoadNorm = Norm(nDim, Load_Dir_Local); - su2double CurrentTime=config->GetCurrent_DynTime(); + su2double CurrentTime=config->GetCurrent_UnstTime(); su2double Ramp_Time = config->GetRamp_Time(); su2double ModAmpl = Compute_LoadCoefficient(CurrentTime, Ramp_Time, config); @@ -2036,7 +2039,7 @@ su2double CFEASolver::Compute_LoadCoefficient(su2double CurrentTime, su2double R /*--- This offset introduces the ramp load in dynamic cases starting from the restart point. ---*/ bool offset = (restart && fsi && (!stat_fsi)); - su2double DeltaT = config->GetDelta_DynTime(); + su2double DeltaT = config->GetDelta_UnstTime(); su2double OffsetTime = offset? DeltaT * (config->GetRestart_Iter()-1) : su2double(0.0); /*--- Polynomial functions from https://en.wikipedia.org/wiki/Smoothstep ---*/ @@ -2507,7 +2510,7 @@ void CFEASolver::Solve_System(CGeometry *geometry, CConfig *config) { void CFEASolver::PredictStruct_Displacement(CGeometry *geometry, const CConfig *config) { const unsigned short predOrder = config->GetPredictorOrder(); - const su2double Delta_t = config->GetDelta_DynTime(); + const su2double Delta_t = config->GetDelta_UnstTime(); const bool dynamic = config->GetTime_Domain(); if(predOrder > 2 && rank == MASTER_NODE) @@ -2575,26 +2578,26 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, const CConfig *c if (RelaxMethod_FSI == BGS_RELAXATION::NONE) { - SetWAitken_Dyn(1.0); + WAitken_Dyn = 1.0; } else if (RelaxMethod_FSI == BGS_RELAXATION::FIXED) { - SetWAitken_Dyn(config->GetAitkenStatRelax()); + WAitken_Dyn = config->GetAitkenStatRelax(); } else if (RelaxMethod_FSI == BGS_RELAXATION::AITKEN) { if (iOuterIter == 0) { - WAitkDyn_tn1 = GetWAitken_Dyn_tn1(); + WAitkDyn_tn1 = WAitken_Dyn_tn1; WAitkDyn_Max = config->GetAitkenDynMaxInit(); WAitkDyn_Min = config->GetAitkenDynMinInit(); WAitkDyn = min(WAitkDyn_tn1, WAitkDyn_Max); WAitkDyn = max(WAitkDyn, WAitkDyn_Min); - SetWAitken_Dyn(WAitkDyn); + WAitken_Dyn = WAitkDyn; } else { @@ -2625,7 +2628,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, const CConfig *c SU2_MPI::Allreduce(&sbuf_numAitk, &rbuf_numAitk, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); SU2_MPI::Allreduce(&sbuf_denAitk, &rbuf_denAitk, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); - WAitkDyn = GetWAitken_Dyn(); + WAitkDyn = WAitken_Dyn; if (rbuf_denAitk > EPS) { WAitkDyn = - 1.0 * WAitkDyn * rbuf_numAitk / rbuf_denAitk ; @@ -2634,7 +2637,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, const CConfig *c WAitkDyn = max(WAitkDyn, 0.1); WAitkDyn = min(WAitkDyn, 1.0); - SetWAitken_Dyn(WAitkDyn); + WAitken_Dyn = WAitkDyn; } @@ -2647,7 +2650,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, const CConfig *c void CFEASolver::SetAitken_Relaxation(CGeometry *geometry, const CConfig *config) { - const su2double WAitken = GetWAitken_Dyn(); + const su2double WAitken = WAitken_Dyn; const bool dynamic = config->GetTime_Domain(); /*--- To nPoint to avoid communication. ---*/ @@ -3107,8 +3110,8 @@ void CFEASolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *c /*--- Detect a wrong solution file. ---*/ if (counter != nPointDomain) { - SU2_MPI::Error(string("The solution file ") + filename + string(" doesn't match with the mesh file!\n") + - string("It could be empty lines at the end of the file."), CURRENT_FUNCTION); + SU2_MPI::Error("The solution file " + filename + " doesn't match with the mesh file!\n" + "It could be empty lines at the end of the file.", CURRENT_FUNCTION); } /*--- MPI. If dynamic, we also need to communicate the old solution. ---*/ @@ -3116,7 +3119,9 @@ void CFEASolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *c InitiateComms(geometry[MESH_0], config, SOLUTION_FEA); CompleteComms(geometry[MESH_0], config, SOLUTION_FEA); - if (dynamic) nodes->Set_Solution_time_n(); + /*--- It's important to not push back the solution when this function is used to load solutions for + * unsteady discrete adjoints, otherwise we overwrite one of the two solutions needed. ---*/ + if (dynamic && val_update_geo) nodes->Set_Solution_time_n(); if (fluid_structure) { for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { diff --git a/SU2_CFD/src/variables/CFEABoundVariable.cpp b/SU2_CFD/src/variables/CFEABoundVariable.cpp index d61f6b2e524c..012498896ad2 100644 --- a/SU2_CFD/src/variables/CFEABoundVariable.cpp +++ b/SU2_CFD/src/variables/CFEABoundVariable.cpp @@ -69,9 +69,10 @@ void CFEABoundVariable::Clear_FlowTraction() { FlowTraction.setConstant(0.0); } void CFEABoundVariable::Clear_SurfaceLoad_Res() { Residual_Ext_Surf.setConstant(0.0); } -void CFEABoundVariable::RegisterFlowTraction() { +void CFEABoundVariable::RegisterFlowTraction(bool reset) { if (!fsi_analysis) return; for (unsigned long iVertex = 0; iVertex < FlowTraction.rows(); iVertex++) for (unsigned long iVar = 0; iVar < nVar; iVar++) - AD::RegisterInput(FlowTraction(iVertex,iVar)); + if (reset) AD::ResetInput(FlowTraction(iVertex,iVar)); + else AD::RegisterInput(FlowTraction(iVertex,iVar)); } diff --git a/SU2_PY/SU2/io/historyMap.py b/SU2_PY/SU2/io/historyMap.py index aca169155e24..a477478359c0 100644 --- a/SU2_PY/SU2/io/historyMap.py +++ b/SU2_PY/SU2/io/historyMap.py @@ -1186,10 +1186,10 @@ "HEADER": "Sens_AoA", "TYPE": "COEFFICIENT", }, - "SENS_E": { + "SENS_E_0": { "DESCRIPTION": "d Objective / d Elasticity modulus", "GROUP": "SENSITIVITY", - "HEADER": "Sens[E]", + "HEADER": "Sens[E_0]", "TYPE": "DEFAULT", }, "SENS_GEO": { @@ -1207,10 +1207,10 @@ "HEADER": "Sens_Mach", "TYPE": "COEFFICIENT", }, - "SENS_NU": { + "SENS_NU_0": { "DESCRIPTION": "d Objective / d Poisson ratio", "GROUP": "SENSITIVITY", - "HEADER": "Sens[Nu]", + "HEADER": "Sens[Nu_0]", "TYPE": "DEFAULT", }, "SENS_PRESS": { diff --git a/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg b/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg index 6b58ea8e0267..a484a26c79b6 100755 --- a/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg +++ b/TestCases/disc_adj_fsi/Airfoil_2d/config.cfg @@ -13,6 +13,6 @@ MESH_FILENAME= mesh.su2 OBJECTIVE_FUNCTION= CUSTOM_OBJFUNC -SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], LINSOL_RESIDUAL[0], SENS_E[1], SENS_NU[1] +SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], AVG_BGS_RES[1], LINSOL_RESIDUAL[0], SENS_E_0[1], SENS_NU_0[1] %WRT_ZONE_CONV=YES diff --git a/TestCases/disc_adj_fsi/config.cfg b/TestCases/disc_adj_fsi/config.cfg index 35cf049125dd..044f409302c1 100644 --- a/TestCases/disc_adj_fsi/config.cfg +++ b/TestCases/disc_adj_fsi/config.cfg @@ -9,6 +9,6 @@ MESH_FILENAME= mesh.su2 OBJECTIVE_FUNCTION = REFERENCE_GEOMETRY -SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], RMS_ADJ_DENSITY[0], SENS_E[1], SENS_NU[1] +SCREEN_OUTPUT= OUTER_ITER, AVG_BGS_RES[0], RMS_ADJ_DENSITY[0], SENS_E_0[1], SENS_NU_0[1] %WRT_ZONE_CONV= YES diff --git a/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref b/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref index e543b661a60b..b2e428bed5a7 100644 --- a/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref +++ b/TestCases/disc_adj_fsi/dyn_fsi/grad_dv.opt.ref @@ -1,9 +1,9 @@ INDEX GRAD -0 -3.461460667601001e-03 -1 -1.841786311588662e-03 -2 -7.915536257748965e-04 -3 -2.739622082729719e-04 -4 -2.734869133461108e-04 -5 -7.881162428890194e-04 -6 -1.828978290516676e-03 -7 -3.427219398258311e-03 +0 -3.461460667601000e-03 +1 -1.841786311588663e-03 +2 -7.915536257748967e-04 +3 -2.739622082729717e-04 +4 -2.734869133461104e-04 +5 -7.881162428890206e-04 +6 -1.828978290516677e-03 +7 -3.427219398258316e-03 diff --git a/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg b/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg index c4181d6339d6..1303f5457be4 100644 --- a/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg +++ b/TestCases/fea_fsi/DynBeam_2d/configBeam_2d.cfg @@ -6,7 +6,7 @@ % Date: 2016.02.01 % % File Version 7.5.1 "Blackbird" % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + SOLVER= ELASTICITY MATH_PROBLEM= DIRECT GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS @@ -18,7 +18,7 @@ RESTART_ITER= 1 ELASTICITY_MODULUS=3E7 POISSON_RATIO=0.3 MATERIAL_DENSITY=7854 -FORMULATION_ELASTICITY_2D = PLANE_STRESS +FORMULATION_ELASTICITY_2D= PLANE_STRESS TIME_DOMAIN=YES TIME_STEP=0.01 MAX_TIME= 0.1 @@ -27,17 +27,17 @@ TIME_ITER=7 TIME_DISCRE_FEA= NEWMARK_IMPLICIT NEWMARK_BETA=0.2601 NEWMARK_GAMMA=0.52 -MARKER_CLAMPED = ( left , right ) -MARKER_PRESSURE= ( lower, 0) -MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0) -LINEAR_SOLVER= FGMRES -LINEAR_SOLVER_PREC= LU_SGS -LINEAR_SOLVER_ERROR= 1E-8 -LINEAR_SOLVER_ITER= 1000 +MARKER_CLAMPED= ( left , right ) +MARKER_PRESSURE= ( lower, 0 ) +MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0 ) +LINEAR_SOLVER= CONJUGATE_GRADIENT +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-6 +LINEAR_SOLVER_ITER= 100 MESH_FORMAT= SU2 TABULAR_FORMAT= CSV VOLUME_FILENAME= beam RESTART_FILENAME= restart_beam.dat SOLUTION_FILENAME= solution_beam.dat OUTPUT_WRT_FREQ= 1 -OUTPUT_FILES= (RESTART_ASCII) +OUTPUT_FILES= RESTART_ASCII, PARAVIEW diff --git a/TestCases/fea_fsi/MixElemsKnowles/config.cfg b/TestCases/fea_fsi/MixElemsKnowles/config.cfg index 94e5328e4b9f..ffa4573286ab 100644 --- a/TestCases/fea_fsi/MixElemsKnowles/config.cfg +++ b/TestCases/fea_fsi/MixElemsKnowles/config.cfg @@ -11,7 +11,6 @@ % Physics SOLVER= ELASTICITY MATH_PROBLEM= DIRECT -DYNAMIC_ANALYSIS= NO % % Optimization OBJECTIVE_FUNCTION= REFERENCE_NODE diff --git a/TestCases/fea_fsi/SquareCyl_Beam/config.cfg b/TestCases/fea_fsi/SquareCyl_Beam/config.cfg index ab740fe42558..d4033e5f4b1b 100644 --- a/TestCases/fea_fsi/SquareCyl_Beam/config.cfg +++ b/TestCases/fea_fsi/SquareCyl_Beam/config.cfg @@ -67,7 +67,7 @@ NONLINEAR_FEM_INT_ITER = 20 % -------------------------- DYNAMIC SIMULATION -------------------------------% % -DYNAMIC_ANALYSIS= YES +TIME_DOMAIN= YES TIME_DISCRE_FEA= NEWMARK_IMPLICIT NEWMARK_BETA=0.2601 NEWMARK_GAMMA=0.52 diff --git a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg index 1aee4275ed54..e3b82a9e6011 100644 --- a/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg +++ b/TestCases/fea_fsi/StatBeam_3d/configBeam_3d.cfg @@ -6,7 +6,7 @@ % Date: 2016.02.01 % % File Version 7.5.1 "Blackbird" % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + SOLVER= ELASTICITY MATH_PROBLEM= DIRECT GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS @@ -15,7 +15,6 @@ MESH_FILENAME= meshBeam_3d.su2 ELASTICITY_MODULUS=3E7 POISSON_RATIO=0.3 MATERIAL_DENSITY=7854 -DYNAMIC_ANALYSIS= NO MARKER_CLAMPED = ( left , right ) MARKER_PRESSURE= ( lower, 0 , symleft, 0, symright, 0) MARKER_LOAD= ( upper, 1, 1000, 0, -1, 0) diff --git a/TestCases/fea_fsi/stat_fsi/configFEA.cfg b/TestCases/fea_fsi/stat_fsi/configFEA.cfg index 3ce2b13322c6..0f4e66ad1f7f 100755 --- a/TestCases/fea_fsi/stat_fsi/configFEA.cfg +++ b/TestCases/fea_fsi/stat_fsi/configFEA.cfg @@ -59,7 +59,7 @@ FORMULATION_ELASTICITY_2D = PLANE_STRAIN POISSON_RATIO=0.4 % -------------------------- DYNAMIC SIMULATION -------------------------------% -DYNAMIC_ANALYSIS= NO +TIME_DOMAIN= NO TIME_DISCRE_FEA= NEWMARK_IMPLICIT % -------------------------- STRUCTURAL SOLVER --------------------------------% diff --git a/TestCases/fea_topology/config.cfg b/TestCases/fea_topology/config.cfg index cf12a17353af..021e878895dd 100644 --- a/TestCases/fea_topology/config.cfg +++ b/TestCases/fea_topology/config.cfg @@ -66,7 +66,6 @@ ITER=1 % Physics SOLVER= ELASTICITY MATH_PROBLEM= DISCRETE_ADJOINT -DYNAMIC_ANALYSIS= NO RESTART_SOL= NO % SOLUTION_FILENAME=solution_structure.dat diff --git a/TestCases/fea_topology/quick_start/settings.cfg b/TestCases/fea_topology/quick_start/settings.cfg index 44d8b3eafaff..8a0edb1041a8 100644 --- a/TestCases/fea_topology/quick_start/settings.cfg +++ b/TestCases/fea_topology/quick_start/settings.cfg @@ -1,7 +1,6 @@ % Physics SOLVER= ELASTICITY MATH_PROBLEM= DIRECT -DYNAMIC_ANALYSIS= NO % ITER= 1 % diff --git a/TestCases/fea_topology/quick_start/settings_compliance.cfg b/TestCases/fea_topology/quick_start/settings_compliance.cfg index d9cf275df75c..57241e0615f9 100644 --- a/TestCases/fea_topology/quick_start/settings_compliance.cfg +++ b/TestCases/fea_topology/quick_start/settings_compliance.cfg @@ -1,7 +1,6 @@ % Physics SOLVER= ELASTICITY MATH_PROBLEM= DISCRETE_ADJOINT -DYNAMIC_ANALYSIS= NO % ITER= 1 % diff --git a/TestCases/fea_topology/quick_start/settings_volfrac.cfg b/TestCases/fea_topology/quick_start/settings_volfrac.cfg index a3ec2ebc3cb0..c7da2fe14f52 100644 --- a/TestCases/fea_topology/quick_start/settings_volfrac.cfg +++ b/TestCases/fea_topology/quick_start/settings_volfrac.cfg @@ -1,7 +1,6 @@ % Physics SOLVER= ELASTICITY MATH_PROBLEM= DISCRETE_ADJOINT -DYNAMIC_ANALYSIS= NO % ITER= 1 % diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 6f69075b1991..95c486de5616 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -1281,6 +1281,15 @@ def main(): pywrapper_aeroelastic.unsteady = True test_list.append(pywrapper_aeroelastic) + # Custom FEA load + pywrapper_custom_fea_load = TestCase('pywrapper_custom_fea_load') + pywrapper_custom_fea_load.cfg_dir = "py_wrapper/custom_load_fea" + pywrapper_custom_fea_load.cfg_file = "config.cfg" + pywrapper_custom_fea_load.test_iter = 13 + pywrapper_custom_fea_load.test_vals = [-7.263559, -4.946814, -14.165142, 34, -6.380144, 320.58] + pywrapper_custom_fea_load.command = TestCase.Command("mpirun -np 2", "python", "run.py") + test_list.append(pywrapper_custom_fea_load) + # FSI, 2d pywrapper_fsi2d = TestCase('pywrapper_fsi2d') pywrapper_fsi2d.cfg_dir = "fea_fsi/WallChannel_2d" diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index a67de6139dbf..61c79202fb80 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -276,7 +276,7 @@ def main(): discadj_fsi2.cfg_dir = "disc_adj_fsi/Airfoil_2d" discadj_fsi2.cfg_file = "config.cfg" discadj_fsi2.test_iter = 8 - discadj_fsi2.test_vals = [-4.349377, 0.128475, -1.303589, 7.5407e-09, 2.3244] + discadj_fsi2.test_vals = [-4.349377, 0.192713, -1.303589, 7.5407e-09, 2.3244] discadj_fsi2.test_vals_aarch64 = [-3.479505, 0.127953, -1.303589, 7.5407e-09, 2.3244] discadj_fsi2.tol = 0.00001 test_list.append(discadj_fsi2) @@ -408,6 +408,19 @@ def main(): test_list.append(pywrapper_FEA_AD_FlowLoad) pass_list.append(pywrapper_FEA_AD_FlowLoad.run_test()) + # FEA unsteady AD Load Sensitivity + pywrapper_Unst_FEA_AD = TestCase('pywrapper_Unst_FEA_AD') + pywrapper_Unst_FEA_AD.cfg_dir = "py_wrapper/custom_load_fea" + pywrapper_Unst_FEA_AD.cfg_file = "config.cfg" + pywrapper_Unst_FEA_AD.test_iter = 100 + pywrapper_Unst_FEA_AD.test_vals = [0.256684, 0.256684, 0.319877, 0.320149, -0.184491, -0.184509] + pywrapper_Unst_FEA_AD.command = TestCase.Command("mpirun -n 2", "python", "run_ad.py") + pywrapper_Unst_FEA_AD.timeout = 1600 + pywrapper_Unst_FEA_AD.tol = 0.00001 + pywrapper_Unst_FEA_AD.new_output = False + test_list.append(pywrapper_Unst_FEA_AD) + pass_list.append(pywrapper_Unst_FEA_AD.run_test()) + # Flow AD Mesh Displacement Sensitivity pywrapper_CFD_AD_MeshDisp = TestCase('pywrapper_CFD_AD_MeshDisp') pywrapper_CFD_AD_MeshDisp.cfg_dir = "py_wrapper/disc_adj_flow/mesh_disp_sens" diff --git a/TestCases/py_wrapper/custom_load_fea/config.cfg b/TestCases/py_wrapper/custom_load_fea/config.cfg new file mode 100644 index 000000000000..c0352190d77d --- /dev/null +++ b/TestCases/py_wrapper/custom_load_fea/config.cfg @@ -0,0 +1,32 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SU2 configuration file % +% Case description: 2D Beam with custom load via Python wrapper % +% File Version 7.5.1 "Blackbird" % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SOLVER= ELASTICITY +GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS +FORMULATION_ELASTICITY_2D= PLANE_STRESS + +MATERIAL_MODEL= NEO_HOOKEAN +ELASTICITY_MODULUS= 10000 +POISSON_RATIO= 0.3 + +MARKER_CLAMPED= ( x_minus ) +MARKER_FLUID_LOAD= ( x_plus, y_minus, y_plus ) + +LINEAR_SOLVER= CONJUGATE_GRADIENT +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-6 +LINEAR_SOLVER_ITER= 100 + +MESH_FORMAT= RECTANGLE +MESH_BOX_SIZE= ( 17, 5, 0 ) +MESH_BOX_LENGTH= ( 0.5, 0.05, 0 ) + +OUTPUT_FILES= RESTART_ASCII, PARAVIEW +SCREEN_OUTPUT= INNER_ITER, RMS_RES, LINSOL, VMS + +INNER_ITER= 20 +CONV_FIELD= REL_RMS_RTOL +CONV_RESIDUAL_MINVAL= -6 diff --git a/TestCases/py_wrapper/custom_load_fea/run.py b/TestCases/py_wrapper/custom_load_fea/run.py new file mode 100755 index 000000000000..24881ef85df8 --- /dev/null +++ b/TestCases/py_wrapper/custom_load_fea/run.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python + +## \file run.py +# \brief FEA case with custom load. +# \version 7.5.1 "Blackbird" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md) +# +# SU2 is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# SU2 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with SU2. If not, see . + +import pysu2 +from mpi4py import MPI + +def main(): + comm = MPI.COMM_WORLD + + # Initialize the corresponding driver of SU2, this includes solver preprocessing. + try: + SU2Driver = pysu2.CSinglezoneDriver('config.cfg', 1, comm) + except TypeError as exception: + print('A TypeError occured in pysu2.CDriver : ', exception) + raise + + # Get the ID of the marker we want to deform. + AllMarkerIDs = SU2Driver.GetMarkerIndices() + MarkerName = 'y_minus' + MarkerID = AllMarkerIDs[MarkerName] if MarkerName in AllMarkerIDs else -1 + + # Number of vertices on the specified marker (per rank). + nVertex = SU2Driver.GetNumberMarkerNodes(MarkerID) if MarkerID >= 0 else 0 + + # Apply a load based on the coordinates. + if nVertex > 0: + MarkerCoords = SU2Driver.MarkerCoordinates(MarkerID) + L = 0.5 + dx = L / 16 # known from mesh settings in this case. + for iVertex in range(nVertex): + x = MarkerCoords(iVertex, 0) + nodalForce = (2 * x / L) * dx + # Half load due to half dx on first and last node. + if abs(x) < 1e-6 or abs(x - L) < 1e-6: + nodalForce = nodalForce / 2 + SU2Driver.SetMarkerCustomFEALoad(MarkerID, iVertex, (0, nodalForce)) + + # Solve. + SU2Driver.StartSolver() + + # Find the tip displacement. + MarkerName = 'x_plus' + MarkerID = AllMarkerIDs[MarkerName] if MarkerName in AllMarkerIDs else -1 + nVertex = SU2Driver.GetNumberMarkerNodes(MarkerID) if MarkerID >= 0 else 0 + Disp = 0 + NodeFound = False + + if nVertex > 0: + MarkerCoords = SU2Driver.MarkerCoordinates(MarkerID) + SolverID = SU2Driver.GetSolverIndices()["FEA"] + Solution = SU2Driver.MarkerSolution(SolverID, MarkerID) + DispID = SU2Driver.GetFEASolutionIndices()["DISPLACEMENT_Y"] + + for iVertex in range(nVertex): + y = MarkerCoords(iVertex, 1) + if abs(y - 0.025) < 1e-6: + Disp = Solution(iVertex, DispID) + NodeFound = True + + if NodeFound: + print(f"Vertical displacement of tip: {Disp}") + # Test the value against expected. + assert abs(Disp / 0.095439 - 1) < 1e-5, "Test FAILED" + + # Finalize the solver and exit cleanly. + SU2Driver.Finalize() + + +if __name__ == '__main__': + main() diff --git a/TestCases/py_wrapper/custom_load_fea/run_ad.py b/TestCases/py_wrapper/custom_load_fea/run_ad.py new file mode 100644 index 000000000000..dd64982b51b4 --- /dev/null +++ b/TestCases/py_wrapper/custom_load_fea/run_ad.py @@ -0,0 +1,250 @@ +#!/usr/bin/env python + +## \file run.py +# \brief Unsteady adjoint FEA case with custom load. +# \version 7.5.1 "Blackbird" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md) +# +# SU2 is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# SU2 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with SU2. If not, see . + +import pysu2ad +from mpi4py import MPI + +common_settings = """ +SOLVER= ELASTICITY +GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS +FORMULATION_ELASTICITY_2D= PLANE_STRESS + +TIME_DOMAIN= YES +TIME_STEP=0.01 +TIME_DISCRE_FEA= NEWMARK_IMPLICIT +NEWMARK_GAMMA= 0.5 +NEWMARK_BETA= 0.25 + +MATERIAL_MODEL= NEO_HOOKEAN +ELASTICITY_MODULUS= 10000 +POISSON_RATIO= 0.3 +MATERIAL_DENSITY= __DENSITY__ + +MARKER_CLAMPED= ( x_minus ) +MARKER_FLUID_LOAD= ( x_plus, y_minus, y_plus ) + +LINEAR_SOLVER= CONJUGATE_GRADIENT +DISCADJ_LIN_SOLVER= CONJUGATE_GRADIENT +LINEAR_SOLVER_PREC= ILU +DISCADJ_LIN_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-5 +LINEAR_SOLVER_ITER= 100 + +MESH_FORMAT= RECTANGLE +MESH_BOX_SIZE= ( 17, 5, 0 ) +MESH_BOX_LENGTH= ( 0.5, __HEIGHT__, 0 ) + +OUTPUT_FILES= RESTART, PARAVIEW +OUTPUT_WRT_FREQ= 1 +OBJECTIVE_FUNCTION= STRESS_PENALTY +STRESS_PENALTY_PARAM= ( 500, 20 ) + +INNER_ITER= 20 +CONV_RESIDUAL_MINVAL= -4 +CONV_STARTITER= 5 +MAX_TIME= 0.2 +TIME_ITER= 21 +""" + +primal_settings = """ +MATH_PROBLEM= DIRECT +SCREEN_OUTPUT= TIME_ITER, CUR_TIME, INNER_ITER, RMS_RES, LINSOL, VMS, STRESS_PENALTY, TAVG_STRESS_PENALTY +HISTORY_OUTPUT= ITER, RMS_RES, STRUCT_COEFF, TAVG_STRUCT_COEFF +CONV_FIELD= REL_RMS_RTOL +""" + +adjoint_settings = """ +MATH_PROBLEM= DISCRETE_ADJOINT +SCREEN_OUTPUT= TIME_ITER, CUR_TIME, INNER_ITER, ADJOINT_DISP_X, ADJOINT_DISP_Y, LINSOL, SENSITIVITY +CONV_FIELD= ADJOINT_DISP_X, ADJOINT_DISP_Y +FEA_ADVANCED_MODE= YES +UNST_ADJOINT_ITER= 21 +ITER_AVERAGE_OBJ= 0 +SOLUTION_FILENAME= restart.dat +VOLUME_OUTPUT= SENSITIVITY_T0 +""" + + +def ApplyLoad(driver, marker_id, peak_load): + """ + Apply a load based on the coordinates and return the derivatives + of the nodal forces with respect to the peak load. + """ + derivatives = [] + if marker_id < 0: return derivatives + + marker_coords = driver.MarkerCoordinates(marker_id) + l = 0.5 + dx = l / 16 # known from mesh settings in this case. + for i_vertex in range(driver.GetNumberMarkerNodes(marker_id)): + x = marker_coords(i_vertex, 0) + nodal_force = (peak_load * x / l) * dx + # Half load due to half dx on first and last node. + if abs(x) < 1e-6 or abs(x - l) < 1e-6: + nodal_force = nodal_force / 2 + driver.SetMarkerCustomFEALoad(marker_id, i_vertex, (0, nodal_force)) + derivatives.append(nodal_force / peak_load) + + return derivatives + + +def RunPrimal(density, peak_load, height): + """ + Runs the primal solver for a given density, peak load, and beam height. + Returns the time average objective function. + """ + comm = MPI.COMM_WORLD + + with open('config_unsteady.cfg', 'w') as f: + f.write(common_settings.replace('__DENSITY__', str(density)).replace('__HEIGHT__', str(height)) + + primal_settings) + + # Initialize the primal driver of SU2, this includes solver preprocessing. + try: + driver = pysu2ad.CSinglezoneDriver('config_unsteady.cfg', 1, comm) + except TypeError as exception: + print('A TypeError occured in pysu2ad.CSinglezoneDriver : ', exception) + raise + + # Get the ID of the marker where the load is applied. + all_marker_ids = driver.GetMarkerIndices() + marker_name = 'y_minus' + marker_id = all_marker_ids[marker_name] if marker_name in all_marker_ids else -1 + + # Apply a load based on the coordinates. + ApplyLoad(driver, marker_id, peak_load) + + # Solve. + driver.StartSolver() + + # Get the time average + tavg_stress_penalty = driver.GetOutputValue('TAVG_STRESS_PENALTY') + + # Finalize the solver and exit cleanly. + driver.Finalize() + + return tavg_stress_penalty + + +def RunAdjoint(density, peak_load, height): + """ + Runs the adjoint solver and returns the sensitivity of the objective function to the peak + load, to the material density, and to the beam height. + """ + comm = MPI.COMM_WORLD + + with open('config_unsteady_ad.cfg', 'w') as f: + f.write(common_settings.replace('__DENSITY__', str(density)).replace('__HEIGHT__', str(height)) + + adjoint_settings) + + # Initialize the adjoint driver of SU2, this includes solver preprocessing. + try: + driver = pysu2ad.CDiscAdjSinglezoneDriver('config_unsteady_ad.cfg', 1, comm) + except TypeError as exception: + print('A TypeError occured in pysu2ad.CDiscAdjSinglezoneDriver : ', exception) + raise + + # Get the ID of the marker where the load is applied. + all_marker_ids = driver.GetMarkerIndices() + marker_name = 'y_minus' + marker_id = all_marker_ids[marker_name] if marker_name in all_marker_ids else -1 + + # Apply the same load that was used in the primal problem. + derivatives = ApplyLoad(driver, marker_id, peak_load) + + n_vertex = driver.GetNumberMarkerNodes(marker_id) if marker_id >= 0 else 0 + load_sens = 0.0 + + # Run the time loop in python to extract sensitivities at each step. + for time_iter in range(driver.GetNumberTimeIter()): + # Preprocess adjoint iteration (AD recording). + driver.Preprocess(time_iter) + + # Run one time iteration. + driver.Run() + driver.Postprocess() + driver.Update() + + # Accumulate load sensitivies (the solver doesn't accumulate + # these for when they are used for FSI adjoints). + for i_vertex in range(n_vertex): + load_sens += derivatives[i_vertex] * driver.GetMarkerFEALoadSensitivity(marker_id, i_vertex)[1] + + # Monitor the solver and output solution to file if required. + driver.Monitor(time_iter) + driver.Output(time_iter) + + rho_sens = driver.GetOutputValue('SENS_RHO_0') + + height_sens = 0.0 + sensitivity = driver.Sensitivity(driver.GetSolverIndices()['ADJ.FEA']) + coords = driver.Coordinates() + + for i_node in range(driver.GetNumberNodes() - driver.GetNumberHaloNodes()): + y = coords(i_node, 1) + dy_dh = y / height + height_sens += dy_dh * sensitivity(i_node, 1) + + # Finalize the solver and exit cleanly. + driver.Finalize() + + return comm.allreduce(load_sens), rho_sens, comm.allreduce(height_sens) + + +def main(): + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + + # Run the primal with 2 loads to compute the sensitivity via finite differences. + obj_pert_height = RunPrimal(1, 2, 0.05001) + obj_pert_load = RunPrimal(1, 2.002, 0.05) + obj_pert_rho = RunPrimal(1.0001, 2, 0.05) + # Run the un-perturbed last to use the restarts for the adjoint. + obj = RunPrimal(1, 2, 0.05) + sens_height_fd = (obj_pert_height - obj) / 0.00001 + sens_load_fd = (obj_pert_load - obj) / 0.002 + sens_rho_fd = (obj_pert_rho - obj) / 0.0001 + + sens_load, sens_rho, sens_height = RunAdjoint(1, 2, 0.05) + + if rank == 0: + print(" Finite Differences\tDiscrete Adjoint") + print(f"Height {sens_height_fd}\t{sens_height}") + print(f"Load {sens_load_fd}\t{sens_load}") + print(f"Rho {sens_rho_fd}\t{sens_rho}") + + assert abs(sens_height / sens_height_fd - 1) < 1e-4, "Error in geometric derivatives." + assert abs(sens_load / sens_load_fd - 1) < 1e-4, "Error in load derivative." + assert abs(sens_rho / sens_rho_fd - 1) < 1e-3, "Error in material derivative." + + # Print results for the regression script to check. + if rank == 0: + print("\n") + print(100, 100, sens_load_fd, sens_load, sens_rho_fd * 10, sens_rho * 10, sens_height_fd / 100, sens_height / 100) + + +if __name__ == '__main__': + main() diff --git a/config_template.cfg b/config_template.cfg index de994a0ce608..a99ec1a89227 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -733,7 +733,7 @@ CONV_NUM_METHOD_SPECIES= SCALAR_UPWIND % Required for 2nd order upwind schemes (NO, YES) MUSCL_SPECIES= NO % -% Slope limiter for species equations (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, BARTH_JESPERSEN, VAN_ALBADA_EDGE) +% Slope limiter for species equations (same as SLOPE_LIMITER_TURB) SLOPE_LIMITER_SPECIES = NONE % % Time discretization for species equations (EULER_IMPLICIT, EULER_EXPLICIT) @@ -1095,36 +1095,37 @@ CUSTOM_OBJFUNC= 'DRAG + 10 * pow(fmax(0.4-LIFT, 0), 2)' % Required for 2nd order upwind schemes (NO, YES) MUSCL_FLOW= YES % -% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, BARTH_JESPERSEN, VAN_ALBADA_EDGE) -% +% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, BARTH_JESPERSEN, VAN_ALBADA_EDGE, +% NISHIKAWA_R3, NISHIKAWA_R4, NISHIKAWA_R5) SLOPE_LIMITER_FLOW= VENKATAKRISHNAN % -% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations. -% Required for 2nd order upwind schemes (NO, YES) +% Same as MUSCL_FLOW but for turbulence. +% MUSCL_TURB= NO % -% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, BARTH_JESPERSEN) +% Slope limiter (same as SLOPE_LIMITER_FLOW except VAN_ALBADA_EDGE) % SLOPE_LIMITER_TURB= VENKATAKRISHNAN % -% Monotonic Upwind Scheme for Conservation Laws (TVD) in the adjoint flow equations. -% Required for 2nd order upwind schemes (NO, YES) +% Same as MUSCL_FLOW but for the continuous adjoint equations. +% MUSCL_ADJFLOW= YES % -% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG, BARTH_JESPERSEN, VAN_ALBADA_EDGE, -% SHARP_EDGES, WALL_DISTANCE, NISHIKAWA_R3, NISHIKAWA_R4, NISHIKAWA_R5) +% Slope limiter (same as SLOPE_LIMITER_FLOW plus SHARP_EDGES, WALL_DISTANCE) +% SLOPE_LIMITER_ADJFLOW= VENKATAKRISHNAN % -% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence adjoint equations. -% Required for 2nd order upwind schemes (NO, YES) +% Same as MUSCL_FLOW but for continuous adjoint turbulence equations. +% MUSCL_ADJTURB= NO % % Slope limiter (see SLOPE_LIMITER_ADJFLOW) +% SLOPE_LIMITER_ADJTURB= VENKATAKRISHNAN % -% Coefficient for the Venkat's limiter (upwind scheme). A larger values decrease -% the extent of limiting, values approaching zero cause -% lower-order approximation to the solution (0.05 by default) +% Coefficient for Venkatakrishnan-type limiters (upwind scheme). +% A larger value decreases the extent of limiting, values approaching zero +% cause lower-order approximation to the solution (0.05 by default) VENKAT_LIMITER_COEFF= 0.05 % % Reference coefficient for detecting sharp edges (3.0 by default).