diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index ab6c1a2f91f6..3719219bfee8 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -992,7 +992,6 @@ class CConfig { su2double RefGeom_Penalty, /*!< \brief Penalty weight value for the reference geometry objective function. */ RefNode_Penalty, /*!< \brief Penalty weight value for the reference node objective function. */ DV_Penalty; /*!< \brief Penalty weight to add a constraint to the total amount of stiffness. */ - bool addCrossTerm; /*!< \brief Evaluates the need to add the cross term when setting the adjoint output. */ unsigned long Nonphys_Points, /*!< \brief Current number of non-physical points in the solution. */ Nonphys_Reconstr; /*!< \brief Current number of non-physical reconstructions for 2nd-order upwinding. */ bool ParMETIS; /*!< \brief Boolean for activating ParMETIS mode (while testing). */ @@ -2085,7 +2084,7 @@ class CConfig { */ unsigned short GetRefGeom_FileFormat(void) const { return RefGeom_FileFormat; } - /*! + /*! * \brief Formulation for 2D elasticity (plane stress - strain) * \return Flag to 2D elasticity model. */ @@ -2097,17 +2096,6 @@ class CConfig { */ bool GetPrestretch(void) const { return Prestretch; } - /*! - * \brief Decide whether it's necessary to add the cross term for adjoint FSI. - * \return TRUE if it's necessary to add the cross term, FALSE otherwise. - */ - bool Add_CrossTerm(void) const { return addCrossTerm; } - - /*! - * \brief Set the boolean addCrossTerm to true or false. - */ - void Set_CrossTerm(bool needCrossTerm) { addCrossTerm = needCrossTerm; } - /*! * \brief Get the name of the file with the element properties for structural problems. * \return Name of the file with the element properties of the structural problem. @@ -2115,9 +2103,9 @@ class CConfig { string GetFEA_FileName(void) const { return FEA_FileName; } /*! - * \brief Determine if advanced features are used from the element-based FEA analysis (experimental feature). - * \return TRUE is experimental, FALSE is the default behaviour. - */ + * \brief Determine if advanced features are used from the element-based FEA analysis (experimental feature). + * \return TRUE is experimental, FALSE is the default behaviour. + */ inline bool GetAdvanced_FEAElementBased(void) const { return FEAAdvancedMode; } /*! @@ -8627,7 +8615,7 @@ class CConfig { * \param[in] val_index - Index corresponding to the load boundary. * \return The pointer to the sine load values. */ - su2double* GetLoad_Sine(void) { return SineLoad_Coeff; } + const su2double* GetLoad_Sine(void) const { return SineLoad_Coeff; } /*! * \brief Get the kind of load transfer method we want to use for dynamic problems diff --git a/Common/include/geometry/elements/CElement.hpp b/Common/include/geometry/elements/CElement.hpp index c3aae53ef255..d8d6b1ee58dd 100644 --- a/Common/include/geometry/elements/CElement.hpp +++ b/Common/include/geometry/elements/CElement.hpp @@ -449,9 +449,10 @@ class CElement { * the latter are needed for compatibility with shape derivatives, there is no problem registering * because inactive variables are ignored. */ - inline void SetPreaccIn_Coords(void) { + inline void SetPreaccIn_Coords(bool nonlinear = true) { AD::SetPreaccIn(RefCoord.data(), nNodes*nDim); - AD::SetPreaccIn(CurrentCoord.data(), nNodes*nDim); + if (nonlinear) + AD::SetPreaccIn(CurrentCoord.data(), nNodes*nDim); } /*! @@ -462,6 +463,13 @@ class CElement { AD::SetPreaccOut(Kt_a.data(), nNodes*nDim); } + /*! + * \brief Register the mass matrix as a pre-accumulation output. + */ + inline void SetPreaccOut_Mab(void) { + AD::SetPreaccOut(Mab.data(), nNodes*nNodes); + } + }; /*! diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 014355c3952c..83f6a66ed98d 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -705,7 +705,13 @@ class CSysMatrix { * \param[in,out] b - The rhs vector (b := b - A_{*,i} * x_i; b_i = x_i). */ template - void EnforceSolutionAtNode(const unsigned long node_i, const OtherType *x_i, CSysVector & b); + void EnforceSolutionAtNode(unsigned long node_i, const OtherType *x_i, CSysVector & b); + + /*! + * \brief Version of EnforceSolutionAtNode for a single degree of freedom. + */ + template + void EnforceSolutionAtDOF(unsigned long node_i, unsigned long iVar, OtherType x_i, CSysVector & b); /*! * \brief Sets the diagonal entries of the matrix as the sum of the blocks in the corresponding column. @@ -718,8 +724,7 @@ class CSysMatrix { * \param[in] alpha - The scaling constant. * \param[in] B - Matrix being. */ - template - void MatrixMatrixAddition(OtherType alpha, const CSysMatrix& B); + void MatrixMatrixAddition(ScalarType alpha, const CSysMatrix& B); /*! * \brief Performs the product of a sparse matrix by a CSysVector. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 4f53145cb13a..725167f7d878 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2020,15 +2020,10 @@ static const MapType DirectDiff_Var_Map = { enum ENUM_RECORDING { - FLOW_CONS_VARS = 1, + SOLUTION_VARIABLES = 1, MESH_COORDS = 2, - COMBINED = 3, - FEA_DISP_VARS = 4, - FLOW_CROSS_TERM = 5, - FEM_CROSS_TERM_GEOMETRY = 6, - GEOMETRY_CROSS_TERM = 7, - ALL_VARIABLES = 8, - MESH_DEFORM = 9 + MESH_DEFORM = 3, + SOLUTION_AND_MESH = 4 }; /*! @@ -2087,7 +2082,6 @@ enum MPI_QUANTITIES { SOLUTION_OLD = 1, /*!< \brief Conservative solution old communication. */ SOLUTION_GRADIENT = 2, /*!< \brief Conservative solution gradient communication. */ SOLUTION_LIMITER = 3, /*!< \brief Conservative solution limiter communication. */ - SOLUTION_DISPONLY = 4, /*!< \brief Solution displacement only communication. */ SOLUTION_PRED = 5, /*!< \brief Solution predicted communication. */ SOLUTION_PRED_OLD = 6, /*!< \brief Solution predicted old communication. */ SOLUTION_GEOMETRY = 7, /*!< \brief Geometry solution communication. */ diff --git a/Common/src/interface_interpolation/CNearestNeighbor.cpp b/Common/src/interface_interpolation/CNearestNeighbor.cpp index 1ab489bd8d94..aeb42bf1135c 100644 --- a/Common/src/interface_interpolation/CNearestNeighbor.cpp +++ b/Common/src/interface_interpolation/CNearestNeighbor.cpp @@ -37,12 +37,7 @@ struct DonorInfo { su2double dist; unsigned pidx; int proc; - DonorInfo(su2double d = 0.0, unsigned i = 0, int p = 0) : - dist(d), pidx(i), proc(p) { } - bool operator< (const DonorInfo& other) const { - /*--- Global index is used as tie-breaker to make sorted order independent of initial. ---*/ - return (dist != other.dist)? (dist < other.dist) : (pidx < other.pidx); - } + DonorInfo(su2double d = 0.0, unsigned i = 0, int p = 0) : dist(d), pidx(i), proc(p) { } }; CNearestNeighbor::CNearestNeighbor(CGeometry ****geometry_container, const CConfig* const* config, unsigned int iZone, @@ -139,8 +134,13 @@ void CNearestNeighbor::SetTransferCoeff(const CConfig* const* config) { } } - /*--- Find k closest points. ---*/ - partial_sort(donorInfo.begin(), donorInfo.begin()+nDonor, donorInfo.end()); + /*--- Find k closest points (need to define the comparator inline or debug build give wrong results). ---*/ + partial_sort(donorInfo.begin(), donorInfo.begin()+nDonor, donorInfo.end(), + [](const DonorInfo& a, const DonorInfo& b) { + /*--- Global index is used as tie-breaker to make sorted order independent of initial. ---*/ + return (a.dist != b.dist)? (a.dist < b.dist) : (a.pidx < b.pidx); + } + ); /*--- Update stats. ---*/ numTarget += 1; diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 018c3c84a870..abc2ea6e69b3 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -1311,6 +1311,45 @@ void CSysMatrix::EnforceSolutionAtNode(const unsigned long node_i, c } +template +template +void CSysMatrix::EnforceSolutionAtDOF(unsigned long node_i, unsigned long iVar, + OtherType x_i, CSysVector & b) { + + for (auto index = row_ptr[node_i]; index < row_ptr[node_i+1]; ++index) { + + const auto node_j = col_ind[index]; + + /*--- Delete row iVar of block j on row i (bij) and ATTEMPT + * to delete column iVar block i on row j (bji). ---*/ + + auto bij = &matrix[index*nVar*nVar]; + auto bji = GetBlock(node_j, node_i); + + /*--- The "attempt" part. ---*/ + if (bji != nullptr) { + for(auto jVar = 0ul; jVar < nVar; ++jVar) { + /*--- Column product. ---*/ + b[node_j*nVar+jVar] -= bji[jVar*nVar+iVar] * x_i; + /*--- Delete entries. ---*/ + bji[jVar*nVar+iVar] = 0.0; + } + } + + /*--- Delete row. ---*/ + for(auto jVar = 0ul; jVar < nVar; ++jVar) + bij[iVar*nVar+jVar] = 0.0; + + /*--- Set the diagonal entry of the block to 1. ---*/ + if (node_j == node_i) + bij[iVar*(nVar+1)] = 1.0; + } + + /*--- Set know solution in rhs vector. ---*/ + b(node_i, iVar) = x_i; + +} + template void CSysMatrix::SetDiagonalAsColumnSum() { @@ -1329,8 +1368,7 @@ void CSysMatrix::SetDiagonalAsColumnSum() { } template -template -void CSysMatrix::MatrixMatrixAddition(OtherType alpha, const CSysMatrix& B) { +void CSysMatrix::MatrixMatrixAddition(ScalarType alpha, const CSysMatrix& B) { /*--- Check that the sparse structure is shared between the two matrices, * comparing pointers is ok as they are obtained from CGeometry. ---*/ @@ -1344,7 +1382,7 @@ void CSysMatrix::MatrixMatrixAddition(OtherType alpha, const CSysMat SU2_OMP_FOR_STAT(omp_light_size) for (auto i = 0ul; i < nnz*nVar*nEqn; ++i) - matrix[i] += PassiveAssign(alpha*B.matrix[i]); + matrix[i] += alpha*B.matrix[i]; } @@ -1403,7 +1441,7 @@ template class CSysMatrix; template void CSysMatrix::InitiateComms(const CSysVector&, CGeometry*, CConfig*, unsigned short) const; template void CSysMatrix::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short) const; template void CSysMatrix::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector&); -template void CSysMatrix::MatrixMatrixAddition(su2double, const CSysMatrix&); +template void CSysMatrix::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double, CSysVector&); #ifdef CODI_REVERSE_TYPE template class CSysMatrix; @@ -1413,6 +1451,6 @@ template void CSysMatrix::CompleteComms(CSysVector::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short) const; template void CSysMatrix::EnforceSolutionAtNode(unsigned long, const passivedouble*, CSysVector&); template void CSysMatrix::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector&); -template void CSysMatrix::MatrixMatrixAddition(passivedouble, const CSysMatrix&); -template void CSysMatrix::MatrixMatrixAddition(su2double, const CSysMatrix&); +template void CSysMatrix::EnforceSolutionAtDOF(unsigned long, unsigned long, passivedouble, CSysVector&); +template void CSysMatrix::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double, CSysVector&); #endif diff --git a/SU2_CFD/include/iteration_structure.hpp b/SU2_CFD/include/iteration_structure.hpp index 5291eb00380f..65a83789e7f8 100644 --- a/SU2_CFD/include/iteration_structure.hpp +++ b/SU2_CFD/include/iteration_structure.hpp @@ -296,12 +296,6 @@ class CIteration { unsigned short iZone, unsigned short iInst){} - virtual void InitializeAdjoint_CrossTerm(CSolver *****solver, - CGeometry ****geometry, - CConfig **config, - unsigned short iZone, - unsigned short iInst){} - virtual void RegisterInput(CSolver *****solver, CGeometry ****geometry, CConfig** config, @@ -1195,7 +1189,7 @@ class CDiscAdjFluidIteration : public CIteration { * \param[in] config - Definition of the particular problem. * \param[in] iZone - Index of the zone. * \param[in] iInst - Index of the instance. - * \param[in] kind_recording - Kind of recording, either FLOW_CONS_VARS or MESH_COORDS + * \param[in] kind_recording - Kind of recording, either FLOW_SOLUTION_VARIABLES or MESH_COORDS */ void RegisterInput(CSolver *****solver, CGeometry ****geometry, @@ -1219,16 +1213,6 @@ class CDiscAdjFluidIteration : public CIteration { unsigned short iZone, unsigned short iInst); - /*! - * \brief Initializes the adjoints of the output variables of the meanflow iteration - without the contribution of the objective function - * \param[in] solver - Container vector with all the solutions. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] iZone - Index of the zone. - * \param[in] iInst - Index of the instance. - */ - void InitializeAdjoint_CrossTerm(CSolver *****solver, CGeometry ****geometry, CConfig **config, unsigned short iZone, unsigned short iInst); - /*! * \brief Record a single iteration of the direct mean flow system. * \param[in] output - Pointer to the COutput class. @@ -1499,16 +1483,6 @@ class CDiscAdjFEAIteration : public CIteration { */ void InitializeAdjoint(CSolver *****solver, CGeometry ****geometry, CConfig** config, unsigned short iZone, unsigned short iInst); - /*! - * \brief Initializes the adjoints of the output variables of the FEM iteration - without the contribution of the objective function - * \param[in] solver - Container vector with all the solutions. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - * \param[in] iZone - Index of the zone. - * \param[in] iInst - Index of the zone. - */ - void InitializeAdjoint_CrossTerm(CSolver *****solver, CGeometry ****geometry, CConfig **config, unsigned short iZone, unsigned short iInst); - /*! * \brief Record a single iteration of the direct FEM system. * \param[in] output - Pointer to the COutput class. diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 02b6f463ae4f..2935d2c96e92 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -37,8 +37,10 @@ */ class CFEASolver : public CSolver { protected: - enum : size_t {MAXNNODE = 8}; + enum : size_t {MAXNNODE_2D = 4}; + enum : size_t {MAXNNODE_3D = 8}; enum : size_t {MAXNVAR = 3}; + enum : size_t {MAXNDIM = 3}; enum : size_t {OMP_MIN_SIZE = 32}; enum : size_t {OMP_MAX_SIZE = 512}; @@ -48,23 +50,24 @@ class CFEASolver : public CSolver { unsigned short *iElem_iDe = nullptr; /*!< \brief For DE cases, ID of the region considered for each iElem. */ - su2double a_dt[9]; /*!< \brief Integration constants. */ + su2double a_dt[9]; /*!< \brief Time integration constants. */ - su2double Conv_Ref[3]; /*!< \brief Reference values for convergence check: DTOL, RTOL, ETOL */ - su2double Conv_Check[3]; /*!< \brief Current values for convergence check: DTOL, RTOL, ETOL */ + su2double Conv_Check[3]; /*!< \brief Current values for convergence check: UTOL, RTOL, ETOL. */ su2double FSI_Conv[2]; /*!< \brief Values to check the convergence of the FSI problem. */ - su2double loadIncrement; /*!< \brief Coefficient that determines the amount of load which is applied */ + unsigned long idxIncrement; /*!< \brief Index of the current load increment */ + su2double loadIncrement; /*!< \brief Coefficient that determines the amount of load which is applied. */ - su2double WAitken_Dyn; /*!< \brief Aitken's dynamic coefficient */ - su2double WAitken_Dyn_tn1; /*!< \brief Aitken's dynamic coefficient in the previous iteration */ + su2double WAitken_Dyn; /*!< \brief Aitken's dynamic coefficient. */ + su2double WAitken_Dyn_tn1; /*!< \brief Aitken's dynamic coefficient in the previous iteration. */ - su2double PenaltyValue; /*!< \brief Penalty value to maintain total stiffness constant */ + su2double PenaltyValue; /*!< \brief Penalty value to maintain total stiffness constant. */ su2double Total_OFRefGeom; /*!< \brief Total Objective Function: Reference Geometry. */ su2double Total_OFRefNode; /*!< \brief Total Objective Function: Reference Node. */ su2double Total_OFVolFrac; /*!< \brief Total Objective Function: Volume fraction (topology optimization). */ su2double Total_OFCompliance; /*!< \brief Total Objective Function: Compliance (topology optimization). */ + su2double Total_OFCombo = 0.0; /*!< \brief One of the above, for output/history purposes. */ su2double Global_OFRefGeom; /*!< \brief Global Objective Function (added over time steps): Reference Geometry. */ su2double Global_OFRefNode; /*!< \brief Global Objective Function (added over time steps): Reference Node. */ @@ -75,6 +78,19 @@ class CFEASolver : public CSolver { su2double RelaxCoeff; /*!< \brief Relaxation coefficient . */ su2double FSI_Residual; /*!< \brief FSI residual. */ + CSysVector TimeRes_Aux; /*!< \brief Auxiliary vector for adding mass and damping contributions to the residual. */ + CSysVector TimeRes; /*!< \brief Vector for adding mass and damping contributions to the residual */ + CSysVector LinSysReact; /*!< \brief Vector to store the residual before applying the BCs */ + +#ifndef CODI_FORWARD_TYPE + CSysMatrix MassMatrix; /*!< \brief Sparse structure for storing the mass matrix. */ +#else + CSysMatrix MassMatrix; +#endif + + CElement*** element_container = nullptr; /*!< \brief Vector which the define the finite element structure for each problem. */ + CProperty** element_properties = nullptr; /*!< \brief Vector which stores the properties of each element */ + #ifdef HAVE_OMP vector > ElemColoring; /*!< \brief Element colors. */ bool LockStrategy = false; /*!< \brief Whether to use an OpenMP lock to guard updates of the Jacobian. */ @@ -163,7 +179,7 @@ class CFEASolver : public CSolver { * \brief Compute constants for time integration. * \param[in] config - Definition of the particular problem. */ - void Compute_IntegrationConstants(CConfig *config); + void Compute_IntegrationConstants(const CConfig *config); /*! * \brief Write the forward mode gradient to file. @@ -174,7 +190,7 @@ class CFEASolver : public CSolver { * \param[in] der - Value of the derivative. * \param[in] der_avg - Time-averaged value of the derivative. */ - void OutputForwardModeGradient(CConfig *config, + void OutputForwardModeGradient(const CConfig *config, bool newFile, su2double fun, su2double fun_avg, @@ -182,19 +198,6 @@ class CFEASolver : public CSolver { su2double der_avg) const; public: - - CSysVector TimeRes_Aux; /*!< \brief Auxiliary vector for adding mass and damping contributions to the residual. */ - CSysVector TimeRes; /*!< \brief Vector for adding mass and damping contributions to the residual */ - CSysVector LinSysReact; /*!< \brief Vector to store the residual before applying the BCs */ - - CSysVector LinSysSol_Adj; /*!< \brief Vector to store the solution of the adjoint problem */ - CSysVector LinSysRes_Adj; /*!< \brief Vector to store the residual of the adjoint problem */ - - CSysMatrix MassMatrix; /*!< \brief Sparse structure for storing the mass matrix. */ - - CElement*** element_container = nullptr; /*!< \brief Vector which the define the finite element structure for each problem. */ - CProperty** element_properties = nullptr; /*!< \brief Vector which stores the properties of each element */ - /*! * \brief Constructor of the class. */ @@ -241,18 +244,6 @@ class CFEASolver : public CSolver { CConfig *config, unsigned long TimeIter) override; - /*! - * \brief Reset the initial condition for the FEM structural problem. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container with all the solutions. - * \param[in] config - Definition of the particular problem. - * \param[in] ExtIter - External iteration. - */ - void ResetInitialCondition(CGeometry **geometry, - CSolver ***solver_container, - CConfig *config, - unsigned long TimeIter) override; - /*! * \brief Compute the time step for solving the FEM equations. * \param[in] geometry - Geometrical definition of the problem. @@ -273,7 +264,7 @@ class CFEASolver : public CSolver { * \param[in] indexNode - Index of the node. * \param[in] iDim - Dimension required. */ - inline virtual su2double Get_ValCoord(CGeometry *geometry, + inline virtual su2double Get_ValCoord(const CGeometry *geometry, unsigned long indexNode, unsigned short iDim) const { return geometry->node[indexNode]->GetCoord(iDim); @@ -282,125 +273,121 @@ class CFEASolver : public CSolver { /*! * \brief Compute the stiffness matrix of the problem. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ void Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, - CConfig *config) final; + const CConfig *config) final; /*! * \brief Compute the stiffness matrix of the problem and the nodal stress terms at the same time. * \note This is more efficient for full Newton Raphson. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ void Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumerics **numerics, - CConfig *config) final; + const CConfig *config) final; /*! * \brief Compute the mass matrix of the problem. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ void Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, - CConfig *config) final; + const CConfig *config) final; /*! * \brief Compute the mass residual of the problem. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ void Compute_MassRes(CGeometry *geometry, CNumerics **numerics, - CConfig *config) final; + const CConfig *config) final; /*! * \brief Compute the nodal stress terms and add them to the residual. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ void Compute_NodalStressRes(CGeometry *geometry, CNumerics **numerics, - CConfig *config) final; + const CConfig *config) final; /*! * \brief Compute the stress at the nodes for output purposes. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ void Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, - CConfig *config) final; + const CConfig *config); /*! * \brief Compute the dead loads. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ void Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, - CConfig *config) final; + const CConfig *config) final; /*! * \brief Clamped boundary conditions. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. */ void BC_Clamped(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final; /*! * \brief Enforce the solution to be 0 in the clamped nodes - Avoids accumulation of numerical error. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. * \param[in] val_marker - Surface marker where the boundary condition is applied. */ void BC_Clamped_Post(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final; /*! - * \brief A virtual member. + * \brief Symmetry boundary conditions. * \param[in] geometry - Geometrical definition of the problem. * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. * \param[in] val_marker - Surface marker where the boundary condition is applied. */ - - void BC_DispDir(CGeometry *geometry, - CNumerics *numerics, - CConfig *config, - unsigned short val_marker) final; + void BC_Sym_Plane(CGeometry *geometry, + CNumerics *numerics, + const CConfig *config, + unsigned short val_marker) final; /*! - * \brief Impose a displacement (constraint) boundary condition. + * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. * \param[in] val_marker - Surface marker where the boundary condition is applied. */ - inline void BC_Normal_Displacement(CGeometry *geometry, - CNumerics *numerics, - CConfig *config, - unsigned short val_marker) final { } + void BC_DispDir(CGeometry *geometry, + CNumerics *numerics, + const CConfig *config, + unsigned short val_marker) final; /*! * \brief Impose a load boundary condition normal to the boundary. @@ -411,7 +398,7 @@ class CFEASolver : public CSolver { */ void BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final; /*! @@ -423,7 +410,7 @@ class CFEASolver : public CSolver { */ void BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final; /*! @@ -435,7 +422,7 @@ class CFEASolver : public CSolver { */ inline void BC_Sine_Load(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final { } /*! @@ -447,7 +434,7 @@ class CFEASolver : public CSolver { */ void BC_Damper(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final; /*! @@ -459,7 +446,7 @@ class CFEASolver : public CSolver { */ void BC_Deforming(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final; /*! @@ -469,85 +456,60 @@ class CFEASolver : public CSolver { */ void Integrate_FSI_Loads(CGeometry *geometry, const CConfig *config); - /*! - * \brief Update the solution using an implicit solver. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] config - Definition of the particular problem. - */ - inline void ImplicitEuler_Iteration(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final { }; - /*! * \brief Iterate using an implicit Newmark solver. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Numerical methods. * \param[in] config - Definition of the particular problem. */ void ImplicitNewmark_Iteration(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + CNumerics **numerics, + const CConfig *config) final; /*! * \brief Update the solution using an implicit Newmark solver. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ - void ImplicitNewmark_Update(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + void ImplicitNewmark_Update(CGeometry *geometry, CConfig *config) final; /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ - void ImplicitNewmark_Relaxation(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + void ImplicitNewmark_Relaxation(CGeometry *geometry, CConfig *config) final; /*! * \brief Iterate using an implicit Generalized Alpha solver. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Numerical methods. * \param[in] config - Definition of the particular problem. */ void GeneralizedAlpha_Iteration(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + CNumerics **numerics, + const CConfig *config) final; /*! * \brief Update the solution using an implicit Generalized Alpha solver. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ - void GeneralizedAlpha_UpdateDisp(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + void GeneralizedAlpha_UpdateDisp(CGeometry *geometry, CConfig *config) final; /*! * \brief Update the solution using an implicit Generalized Alpha solver. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ - void GeneralizedAlpha_UpdateSolution(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + void GeneralizedAlpha_UpdateSolution(CGeometry *geometry, CConfig *config) final; /*! * \brief Update the solution using an implicit Generalized Alpha solver. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ - void GeneralizedAlpha_UpdateLoads(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + void GeneralizedAlpha_UpdateLoads(CGeometry *geometry, const CConfig *config) final; /*! * \brief Postprocessing. @@ -607,6 +569,12 @@ class CFEASolver : public CSolver { */ inline su2double GetTotal_OFCompliance(void) const final { return Total_OFCompliance; } + /*! + * \brief Retrieve the value of the combined objective function + * \note For now there is no combination, this is just a seletion. + */ + inline su2double GetTotal_ComboObj(void) const final { return Total_OFCombo; } + /*! * \brief Determines whether there is an element-based file or not. * \return Bool that defines whether the solution has an element-based file or not @@ -672,85 +640,55 @@ class CFEASolver : public CSolver { /*! * \brief Predictor for structural displacements based on previous iterations - * \param[in] fea_geometry - Geometrical definition of the problem. - * \param[in] fea_grid_movement - Geometrical definition of the problem. - * \param[in] fea_config - Geometrical definition of the problem. - * \param[in] flow_geometry - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Configuration of the problem. */ - void PredictStruct_Displacement(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) final; + void PredictStruct_Displacement(CGeometry *geometry, CConfig *config) final; /*! * \brief Computation of Aitken's coefficient. - * \param[in] fea_geometry - Geometrical definition of the problem. - * \param[in] fea_config - Geometrical definition of the problem. - * \param[in] fea_geometry - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] iOuterIter - Current outer iteration. */ - void ComputeAitken_Coefficient(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution, + void ComputeAitken_Coefficient(CGeometry *geometry, + CConfig *config, unsigned long iOuterIter) final; /*! * \brief Aitken's relaxation of the solution. - * \param[in] fea_geometry - Geometrical definition of the problem. - * \param[in] fea_config - Geometrical definition of the problem. - * \param[in] fea_geometry - Definition of the particular problem. - */ - void SetAitken_Relaxation(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) final; - - /*! - * \brief Aitken's relaxation of the solution. - * \param[in] fea_geometry - Geometrical definition of the problem. - * \param[in] fea_config - Geometrical definition of the problem. - * \param[in] fea_geometry - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. */ - void Update_StructSolution(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) final; + void SetAitken_Relaxation(CGeometry *geometry, CConfig *config) final; /*! * \brief Compute the objective function for a reference geometry * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ - void Compute_OFRefGeom(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + void Compute_OFRefGeom(CGeometry *geometry, const CConfig *config) final; /*! * \brief Compute the objective function for a reference node * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ - void Compute_OFRefNode(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + void Compute_OFRefNode(CGeometry *geometry, const CConfig *config) final; /*! * \brief Compute the objective function for a volume fraction * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ - void Compute_OFVolFrac(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + void Compute_OFVolFrac(CGeometry *geometry, const CConfig *config) final; /*! * \brief Compute the compliance objective function * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ - void Compute_OFCompliance(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) final; + void Compute_OFCompliance(CGeometry *geometry, const CConfig *config) final; /*! * \brief Compute the penalty due to the stiffness increase @@ -762,7 +700,7 @@ class CFEASolver : public CSolver { void Stiffness_Penalty(CGeometry *geometry, CSolver **solver_container, CNumerics **numerics_container, - CConfig *config) final; + CConfig *config); /*! * \brief Get the value of the FSI convergence. @@ -805,7 +743,10 @@ class CFEASolver : public CSolver { * \brief Set the value of the load increment for nonlinear structural analysis * \param[in] Value of the coefficient */ - inline void SetLoad_Increment(su2double val_loadIncrement) final { loadIncrement = val_loadIncrement; } + inline void SetLoad_Increment(unsigned long iInc, su2double loadInc) final { + idxIncrement = iInc; + loadIncrement = loadInc; + } /*! * \brief Get the value of the load increment for nonlinear structural analysis @@ -841,7 +782,7 @@ class CFEASolver : public CSolver { */ su2double Compute_LoadCoefficient(su2double CurrentTime, su2double RampTime, - CConfig *config) final; + const CConfig *config) final; /*! * \brief A virtual member. @@ -865,6 +806,6 @@ class CFEASolver : public CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void FilterElementDensities(CGeometry *geometry, CConfig *config); + void FilterElementDensities(CGeometry *geometry, const CConfig *config); }; diff --git a/SU2_CFD/include/solvers/CMeshSolver.hpp b/SU2_CFD/include/solvers/CMeshSolver.hpp index 3dcd474c7e5a..d661cc8acab5 100644 --- a/SU2_CFD/include/solvers/CMeshSolver.hpp +++ b/SU2_CFD/include/solvers/CMeshSolver.hpp @@ -130,7 +130,7 @@ class CMeshSolver final : public CFEASolver { * \param[in] indexNode - Index of the node. * \param[in] iDim - Dimension required. */ - inline su2double Get_ValCoord(CGeometry*, + inline su2double Get_ValCoord(const CGeometry*, unsigned long indexNode, unsigned short iDim) const override { return nodes->GetMesh_Coord(indexNode,iDim); diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index d2bf6e363e6a..bc5b71a6ffbb 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -164,7 +164,6 @@ class CSolver { CSysVector LinSysSol; /*!< \brief vector to store iterative solution of implicit linear system. */ CSysVector LinSysRes; /*!< \brief vector to store iterative residual of implicit linear system. */ - CSysVector LinSysAux; /*!< \brief vector to store iterative residual of implicit linear system. */ #ifndef CODI_FORWARD_TYPE CSysMatrix Jacobian; /*!< \brief Complete sparse Jacobian structure for implicit computations. */ CSysSolve System; /*!< \brief Linear solver/smoother. */ @@ -173,8 +172,6 @@ class CSolver { CSysSolve System; #endif - CSysMatrix StiffMatrix; /*!< \brief Sparse structure for storing the stiffness matrix in Galerkin computations, and grid movement. */ - CSysVector OutputVariables; /*!< \brief vector to store the extra variables to be written. */ string* OutputHeadingNames; /*!< \brief vector of strings to store the headings for the exra variables */ @@ -860,23 +857,21 @@ class CSolver { * \param[in] config - Definition of the particular problem. * \param[in] val_marker - Surface marker where the boundary condition is applied. */ - inline virtual void BC_Clamped(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. * \param[in] val_marker - Surface marker where the boundary condition is applied. */ - inline virtual void BC_Clamped_Post(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! @@ -886,10 +881,21 @@ class CSolver { * \param[in] config - Definition of the particular problem. * \param[in] val_marker - Surface marker where the boundary condition is applied. */ + inline virtual void BC_Sym_Plane(CGeometry *geometry, + CNumerics *numerics, + const CConfig *config, + unsigned short val_marker) { } + /*! + * \brief A virtual member. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] numerics - Description of the numerical method. + * \param[in] config - Definition of the particular problem. + * \param[in] val_marker - Surface marker where the boundary condition is applied. + */ inline virtual void BC_DispDir(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! @@ -899,13 +905,11 @@ class CSolver { * \param[in] config - Definition of the particular problem. * \param[in] val_marker - Surface marker where the boundary condition is applied. */ - inline virtual void BC_Normal_Displacement(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } - /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. @@ -915,7 +919,7 @@ class CSolver { */ inline virtual void BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! @@ -925,10 +929,9 @@ class CSolver { * \param[in] config - Definition of the particular problem. * \param[in] val_marker - Surface marker where the boundary condition is applied. */ - inline virtual void BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! @@ -941,7 +944,7 @@ class CSolver { inline virtual void BC_Sine_Load(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! @@ -953,7 +956,7 @@ class CSolver { */ inline virtual void BC_Damper(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! @@ -963,10 +966,9 @@ class CSolver { * \param[in] config - Definition of the particular problem. * \param[in] val_marker - Surface marker where the boundary condition is applied. */ - inline virtual void BC_Deforming(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! @@ -1522,12 +1524,12 @@ class CSolver { /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Numerical methods. * \param[in] config - Definition of the particular problem. */ inline virtual void ImplicitNewmark_Iteration(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) { } + CNumerics **numerics, + const CConfig *config) { } /*! * \brief A virtual member. @@ -1536,7 +1538,6 @@ class CSolver { * \param[in] config - Definition of the particular problem. */ inline virtual void ImplicitNewmark_Update(CGeometry *geometry, - CSolver **solver_container, CConfig *config) { } /*! @@ -1546,18 +1547,17 @@ class CSolver { * \param[in] config - Definition of the particular problem. */ inline virtual void ImplicitNewmark_Relaxation(CGeometry *geometry, - CSolver **solver_container, CConfig *config) { } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics - Numerical methods. * \param[in] config - Definition of the particular problem. */ inline virtual void GeneralizedAlpha_Iteration(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) { } + CNumerics **numerics, + const CConfig *config) { } /*! * \brief A virtual member. @@ -1566,7 +1566,6 @@ class CSolver { * \param[in] config - Definition of the particular problem. */ inline virtual void GeneralizedAlpha_UpdateDisp(CGeometry *geometry, - CSolver **solver_container, CConfig *config) { } /*! @@ -1576,7 +1575,6 @@ class CSolver { * \param[in] config - Definition of the particular problem. */ inline virtual void GeneralizedAlpha_UpdateSolution(CGeometry *geometry, - CSolver **solver_container, CConfig *config) { } /*! @@ -1586,8 +1584,7 @@ class CSolver { * \param[in] config - Definition of the particular problem. */ inline virtual void GeneralizedAlpha_UpdateLoads(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. @@ -3580,56 +3577,28 @@ class CSolver { /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container with all the solutions. - * \param[in] config - Definition of the particular problem. - * \param[in] ExtIter - External iteration. - */ - inline virtual void ResetInitialCondition(CGeometry **geometry, - CSolver ***solver_container, - CConfig *config, - unsigned long TimeIter) { } - - /*! - * \brief A virtual member. - * \param[in] fea_geometry - Geometrical definition of the problem. - * \param[in] fea_config - Geometrical definition of the problem. - * \param[in] fea_geometry - Definition of the particular problem. + * \param[in] config - Definition of the problem. */ - inline virtual void PredictStruct_Displacement(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) { } + inline virtual void PredictStruct_Displacement(CGeometry *geometry, + CConfig *config) { } /*! * \brief A virtual member. - * \param[in] fea_geometry - Geometrical definition of the problem. - * \param[in] fea_config - Geometrical definition of the problem. - * \param[in] fea_geometry - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the problem. + * \param[in] iOuterIter - Current outer iteration. */ - inline virtual void ComputeAitken_Coefficient(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution, + inline virtual void ComputeAitken_Coefficient(CGeometry *geometry, + CConfig *config, unsigned long iOuterIter) { } - - /*! - * \brief A virtual member. - * \param[in] fea_geometry - Geometrical definition of the problem. - * \param[in] fea_config - Geometrical definition of the problem. - * \param[in] fea_geometry - Definition of the particular problem. - */ - inline virtual void SetAitken_Relaxation(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) { } - /*! * \brief A virtual member. - * \param[in] fea_geometry - Geometrical definition of the problem. - * \param[in] fea_config - Geometrical definition of the problem. - * \param[in] fea_geometry - Definition of the particular problem. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. */ - inline virtual void Update_StructSolution(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) { } + inline virtual void SetAitken_Relaxation(CGeometry *geometry, + CConfig *config) { } /*! * \brief A virtual member. @@ -3696,54 +3665,34 @@ class CSolver { /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ inline virtual void Compute_OFRefGeom(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ inline virtual void Compute_OFRefNode(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ inline virtual void Compute_OFVolFrac(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. * \param[in] config - Definition of the particular problem. */ inline virtual void Compute_OFCompliance(CGeometry *geometry, - CSolver **solver_container, - CConfig *config) { } - - /*! - * \brief A virtual member. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] numerics - Description of the numerical method. - * \param[in] config - Definition of the particular problem. - */ - inline virtual void Stiffness_Penalty(CGeometry *geometry, - CSolver **solver_container, - CNumerics **numerics_container, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. @@ -3951,79 +3900,67 @@ class CSolver { */ inline virtual su2double Compute_LoadCoefficient(su2double CurrentTime, su2double RampTime, - CConfig *config) { return 0.0; } + const CConfig *config) { return 0.0; } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ inline virtual void Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ inline virtual void Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumerics **numerics, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ inline virtual void Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ inline virtual void Compute_MassRes(CGeometry *geometry, CNumerics **numerics, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ inline virtual void Compute_NodalStressRes(CGeometry *geometry, CNumerics **numerics, - CConfig *config) { } - - /*! - * \brief A virtual member. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver - Description of the numerical method. - * \param[in] config - Definition of the particular problem. - */ - - inline virtual void Compute_NodalStress(CGeometry *geometry, - CNumerics **numerics, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. - * \param[in] solver_container - Container vector with all the solutions. - * \param[in] solver - Description of the numerical method. + * \param[in] numerics - Description of the numerical method. * \param[in] config - Definition of the particular problem. */ inline virtual void Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. Set the volumetric heat source @@ -4067,15 +4004,16 @@ class CSolver { /*! * \brief A virtual member. - * \param[in] Value of the load increment for nonlinear structural analysis + * \param[in] Index of the increment. + * \param[in] Value of the load increment for nonlinear structural analysis. */ - inline virtual void SetLoad_Increment(su2double val_loadIncrement) { } + inline virtual void SetLoad_Increment(unsigned long iInc, su2double loadInc) { } /*! * \brief A virtual member. * \param[in] Value of the load increment for nonlinear structural analysis */ - inline virtual su2double GetLoad_Increment() const { return 0; } + inline virtual su2double GetLoad_Increment() const { return 0.0; } /*! * \brief A virtual member. diff --git a/SU2_CFD/include/variables/CFEAVariable.hpp b/SU2_CFD/include/variables/CFEAVariable.hpp index 16a083dc7c1d..9611631fb927 100644 --- a/SU2_CFD/include/variables/CFEAVariable.hpp +++ b/SU2_CFD/include/variables/CFEAVariable.hpp @@ -6,7 +6,7 @@ * * SU2 Project Website: https://su2code.github.io * - * The SU2 Project is maintained by the SU2 Foundation + * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md) @@ -337,7 +337,7 @@ class CFEAVariable : public CVariable { /*! * \brief A virtual member. */ - inline su2double *GetPrestretch(unsigned long iPoint) final { return Prestretch[iPoint]; } + inline const su2double *GetPrestretch(unsigned long iPoint) const final { return Prestretch[iPoint]; } /*! * \brief A virtual member. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 5f4850025498..f7289f02299a 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2384,7 +2384,7 @@ class CVariable { /*! * \brief A virtual member. */ - inline virtual su2double *GetPrestretch(unsigned long iPoint) {return nullptr; } + inline virtual const su2double *GetPrestretch(unsigned long iPoint) const {return nullptr; } /*! * \brief A virtual member. diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index b5be9b5a904d..1403f674602e 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -179,7 +179,7 @@ void CDiscAdjMultizoneDriver::Run() { /*--- Evaluate the objective function gradient w.r.t. the solutions of all zones. ---*/ SetRecording(NONE, Kind_Tape::FULL_TAPE, ZONE_0); - SetRecording(FLOW_CONS_VARS, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0); + SetRecording(SOLUTION_VARIABLES, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0); RecordingState = NONE; AD::ClearAdjoints(); @@ -212,9 +212,9 @@ void CDiscAdjMultizoneDriver::Run() { * * To set the tape appropriately, the following recording methods are provided: * (1) NONE: All information from a previous recording is removed. - * (2) FLOW_CONS_VARS: State variables of all solvers in a zone as input. - * (3) MESH_COORDS: Mesh coordinates as input. - * (4) COMBINED: Mesh coordinates and state variables as input. + * (2) SOLUTION_VARIABLES: State variables of all solvers in a zone as input. + * (3) MESH_COORDS / MESH_DEFORM: Mesh coordinates as input. + * (4) SOLUTION_AND_MESH: Mesh coordinates and state variables as input. * * By default, all (state and mesh coordinate variables) will be declared as output, * since it does not change the computational effort. ---*/ @@ -223,9 +223,9 @@ void CDiscAdjMultizoneDriver::Run() { /*--- If we want to set up zone-specific tapes (retape), we do not need to record * here. Otherwise, the whole tape of a coupled run will be created. ---*/ - if (!retape && (RecordingState != FLOW_CONS_VARS)) { + if (!retape && (RecordingState != SOLUTION_VARIABLES)) { SetRecording(NONE, Kind_Tape::FULL_TAPE, ZONE_0); - SetRecording(FLOW_CONS_VARS, Kind_Tape::FULL_TAPE, ZONE_0); + SetRecording(SOLUTION_VARIABLES, Kind_Tape::FULL_TAPE, ZONE_0); } /*-- Start loop over zones. ---*/ @@ -234,7 +234,7 @@ void CDiscAdjMultizoneDriver::Run() { if (retape) { SetRecording(NONE, Kind_Tape::FULL_TAPE, ZONE_0); - SetRecording(FLOW_CONS_VARS, Kind_Tape::ZONE_SPECIFIC_TAPE, iZone); + SetRecording(SOLUTION_VARIABLES, Kind_Tape::ZONE_SPECIFIC_TAPE, iZone); } /*--- Start inner iterations from where we stopped in previous outer iteration. ---*/ @@ -428,9 +428,9 @@ void CDiscAdjMultizoneDriver::SetRecording(unsigned short kind_recording, Kind_T if (rank == MASTER_NODE) { cout << "\n-------------------------------------------------------------------------\n"; switch(kind_recording) { - case NONE: cout << "Clearing the computational graph." << endl; break; - case MESH_COORDS: cout << "Storing computational graph wrt MESH COORDINATES." << endl; break; - case FLOW_CONS_VARS: cout << "Storing computational graph wrt CONSERVATIVE VARIABLES." << endl; break; + case NONE: cout << "Clearing the computational graph." << endl; break; + case MESH_COORDS: cout << "Storing computational graph wrt MESH COORDINATES." << endl; break; + case SOLUTION_VARIABLES: cout << "Storing computational graph wrt CONSERVATIVE VARIABLES." << endl; break; } } @@ -533,7 +533,7 @@ void CDiscAdjMultizoneDriver::DirectIteration(unsigned short iZone, unsigned sho /*--- Print residuals in the first iteration ---*/ - if (rank == MASTER_NODE && kind_recording == FLOW_CONS_VARS) { + if (rank == MASTER_NODE && kind_recording == SOLUTION_VARIABLES) { auto solvers = solver_container[iZone][INST_0][MESH_0]; @@ -734,21 +734,21 @@ void CDiscAdjMultizoneDriver::SetObjFunction(unsigned short kind_recording) { switch(config->GetKind_ObjFunc()) { case REFERENCE_NODE: - solvers[FEA_SOL]->Compute_OFRefNode(geometry, solvers, config); + solvers[FEA_SOL]->Compute_OFRefNode(geometry, config); ObjFunc += solvers[FEA_SOL]->GetTotal_OFRefNode()*Weight_ObjFunc; break; case REFERENCE_GEOMETRY: - solvers[FEA_SOL]->Compute_OFRefGeom(geometry, solvers, config); + solvers[FEA_SOL]->Compute_OFRefGeom(geometry, config); ObjFunc += solvers[FEA_SOL]->GetTotal_OFRefGeom()*Weight_ObjFunc; break; case TOPOL_COMPLIANCE: static_cast(solvers[FEA_SOL])->Integrate_FSI_Loads(geometry, config); - solvers[FEA_SOL]->Compute_OFCompliance(geometry, solvers, config); + solvers[FEA_SOL]->Compute_OFCompliance(geometry, config); ObjFunc += solvers[FEA_SOL]->GetTotal_OFCompliance()*Weight_ObjFunc; break; case VOLUME_FRACTION: case TOPOL_DISCRETENESS: - solvers[FEA_SOL]->Compute_OFVolFrac(geometry, solvers, config); + solvers[FEA_SOL]->Compute_OFVolFrac(geometry, config); ObjFunc += solvers[FEA_SOL]->GetTotal_OFVolFrac()*Weight_ObjFunc; break; @@ -762,14 +762,14 @@ void CDiscAdjMultizoneDriver::SetObjFunction(unsigned short kind_recording) { break; } - if (ObjectiveNotCovered && (rank == MASTER_NODE) && (kind_recording == FLOW_CONS_VARS)) + if (ObjectiveNotCovered && (rank == MASTER_NODE) && (kind_recording == SOLUTION_VARIABLES)) cout << " Objective function not covered in Zone " << iZone << endl; } if (rank == MASTER_NODE) { AD::RegisterOutput(ObjFunc); AD::SetIndex(ObjFunc_Index, ObjFunc); - if (kind_recording == FLOW_CONS_VARS) { + if (kind_recording == SOLUTION_VARIABLES) { cout << " Objective function : " << ObjFunc; if (driver_config->GetWrt_AD_Statistics()){ cout << " (" << ObjFunc_Index << ")\n"; diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp index 18d7da37d96d..3e068ebd4726 100644 --- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp @@ -6,7 +6,7 @@ * * SU2 Project Website: https://su2code.github.io * - * The SU2 Project is maintained by the SU2 Foundation + * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md) @@ -75,7 +75,7 @@ CDiscAdjSinglezoneDriver::CDiscAdjSinglezoneDriver(char* confFile, else direct_iteration = new CFluidIteration(config); if (compressible) direct_output = COutputFactory::createOutput(EULER, config, nDim); else direct_output = COutputFactory::createOutput(INC_EULER, config, nDim); - MainVariables = FLOW_CONS_VARS; + MainVariables = SOLUTION_VARIABLES; if (mesh_def) SecondaryVariables = MESH_DEFORM; else SecondaryVariables = MESH_COORDS; break; @@ -85,7 +85,7 @@ CDiscAdjSinglezoneDriver::CDiscAdjSinglezoneDriver(char* confFile, cout << "Direct iteration: Euler/Navier-Stokes/RANS equation." << endl; direct_iteration = new CFEMFluidIteration(config); direct_output = COutputFactory::createOutput(FEM_EULER, config, nDim); - MainVariables = FLOW_CONS_VARS; + MainVariables = SOLUTION_VARIABLES; SecondaryVariables = MESH_COORDS; break; @@ -94,7 +94,7 @@ CDiscAdjSinglezoneDriver::CDiscAdjSinglezoneDriver(char* confFile, cout << "Direct iteration: elasticity equation." << endl; direct_iteration = new CFEAIteration(config); direct_output = COutputFactory::createOutput(FEM_ELASTICITY, config, nDim); - MainVariables = FEA_DISP_VARS; + MainVariables = SOLUTION_VARIABLES; SecondaryVariables = MESH_COORDS; break; @@ -103,7 +103,7 @@ CDiscAdjSinglezoneDriver::CDiscAdjSinglezoneDriver(char* confFile, cout << "Direct iteration: heat equation." << endl; direct_iteration = new CHeatIteration(config); direct_output = COutputFactory::createOutput(HEAT_EQUATION, config, nDim); - MainVariables = FLOW_CONS_VARS; + MainVariables = SOLUTION_VARIABLES; SecondaryVariables = MESH_COORDS; break; diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index ec76280e1d38..d438da173421 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -254,13 +254,10 @@ void CMultizoneDriver::Preprocess(unsigned long TimeIter) { /*--- Set the initial condition for EULER/N-S/RANS ---------------------------------------------*/ /*--- For FSI, this is set after the mesh has been moved. --------------------------------------*/ - if ((config_container[iZone]->GetKind_Solver() == EULER) || - (config_container[iZone]->GetKind_Solver() == NAVIER_STOKES) || - (config_container[iZone]->GetKind_Solver() == RANS) || - (config_container[iZone]->GetKind_Solver() == INC_EULER) || - (config_container[iZone]->GetKind_Solver() == INC_NAVIER_STOKES) || - (config_container[iZone]->GetKind_Solver() == INC_RANS) ) { - if(!fsi) solver_container[iZone][INST_0][MESH_0][FLOW_SOL]->SetInitialCondition(geometry_container[iZone][INST_0], solver_container[iZone][INST_0], config_container[iZone], TimeIter); + if (!fsi && !config_container[iZone]->GetDiscrete_Adjoint() && config_container[iZone]->GetFluidProblem()) { + solver_container[iZone][INST_0][MESH_0][FLOW_SOL]->SetInitialCondition(geometry_container[iZone][INST_0], + solver_container[iZone][INST_0], + config_container[iZone], TimeIter); } } @@ -270,10 +267,11 @@ void CMultizoneDriver::Preprocess(unsigned long TimeIter) { #endif /*--- Run a predictor step ---*/ - for (iZone = 0; iZone < nZone; iZone++){ + for (iZone = 0; iZone < nZone; iZone++) { if (config_container[iZone]->GetPredictor()) - iteration_container[iZone][INST_0]->Predictor(output_container[iZone], integration_container, geometry_container, solver_container, - numerics_container, config_container, surface_movement, grid_movement, FFDBox, iZone, INST_0); + iteration_container[iZone][INST_0]->Predictor(output_container[iZone], integration_container, geometry_container, + solver_container, numerics_container, config_container, surface_movement, + grid_movement, FFDBox, iZone, INST_0); } /*--- Perform a dynamic mesh update if required. ---*/ diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index d37de0475dab..5b72073a9546 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -331,12 +331,11 @@ void CIntegration::SetStructural_Solver(CGeometry *geometry, CSolver **solver_co case (CD_EXPLICIT): break; case (NEWMARK_IMPLICIT): - if (fsi) solver_container[FEA_SOL]->ImplicitNewmark_Relaxation(geometry, solver_container, config); + if (fsi) solver_container[FEA_SOL]->ImplicitNewmark_Relaxation(geometry, config); break; case (GENERALIZED_ALPHA): - //if (fsi) solver_container[FEA_SOL]->Update_StructSolution(geometry, solver_container, config); - solver_container[FEA_SOL]->GeneralizedAlpha_UpdateSolution(geometry, solver_container, config); - solver_container[FEA_SOL]->GeneralizedAlpha_UpdateLoads(geometry, solver_container, config); + solver_container[FEA_SOL]->GeneralizedAlpha_UpdateSolution(geometry, config); + solver_container[FEA_SOL]->GeneralizedAlpha_UpdateLoads(geometry, config); break; } diff --git a/SU2_CFD/src/integration/CStructuralIntegration.cpp b/SU2_CFD/src/integration/CStructuralIntegration.cpp index ce9b9bd49de8..b73577dba315 100644 --- a/SU2_CFD/src/integration/CStructuralIntegration.cpp +++ b/SU2_CFD/src/integration/CStructuralIntegration.cpp @@ -137,13 +137,13 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * switch (config->GetKind_TimeIntScheme_FEA()) { case (CD_EXPLICIT): - solver_container[MainSolver]->ImplicitNewmark_Iteration(geometry, solver_container, config); + solver_container[MainSolver]->ImplicitNewmark_Iteration(geometry, numerics, config); break; case (NEWMARK_IMPLICIT): - solver_container[MainSolver]->ImplicitNewmark_Iteration(geometry, solver_container, config); + solver_container[MainSolver]->ImplicitNewmark_Iteration(geometry, numerics, config); break; case (GENERALIZED_ALPHA): - solver_container[MainSolver]->GeneralizedAlpha_Iteration(geometry, solver_container, config); + solver_container[MainSolver]->GeneralizedAlpha_Iteration(geometry, numerics, config); break; } @@ -154,16 +154,16 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * case CLAMPED_BOUNDARY: solver_container[MainSolver]->BC_Clamped(geometry, numerics[FEA_TERM], config, iMarker); break; + case SYMMETRY_PLANE: + solver_container[MainSolver]->BC_Sym_Plane(geometry, numerics[FEA_TERM], config, iMarker); + break; case DISP_DIR_BOUNDARY: solver_container[MainSolver]->BC_DispDir(geometry, numerics[FEA_TERM], config, iMarker); break; - case DISPLACEMENT_BOUNDARY: - solver_container[MainSolver]->BC_Normal_Displacement(geometry, numerics[CONV_BOUND_TERM], config, iMarker); - break; } } - /*--- Solver linearized system ---*/ + /*--- Solver linear system ---*/ solver_container[MainSolver]->Solve_System(geometry, config); @@ -171,13 +171,13 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * switch (config->GetKind_TimeIntScheme_FEA()) { case (CD_EXPLICIT): - solver_container[MainSolver]->ImplicitNewmark_Update(geometry, solver_container, config); + solver_container[MainSolver]->ImplicitNewmark_Update(geometry, config); break; case (NEWMARK_IMPLICIT): - solver_container[MainSolver]->ImplicitNewmark_Update(geometry, solver_container, config); + solver_container[MainSolver]->ImplicitNewmark_Update(geometry, config); break; case (GENERALIZED_ALPHA): - solver_container[MainSolver]->GeneralizedAlpha_UpdateDisp(geometry, solver_container, config); + solver_container[MainSolver]->GeneralizedAlpha_UpdateDisp(geometry, config); break; } @@ -191,9 +191,4 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * } } - /*--- Perform the MPI communication of the solution ---*/ - - solver_container[MainSolver]->InitiateComms(geometry, config, SOLUTION_FEA); - solver_container[MainSolver]->CompleteComms(geometry, config, SOLUTION_FEA); - } diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index 972cd245e213..54efa417df01 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -1212,17 +1212,18 @@ void CFEAIteration::Iterate(COutput *output, unsigned short val_iZone, unsigned short val_iInst) { - su2double loadIncrement; + bool StopCalc = false; unsigned long IntIter = 0; - unsigned long TimeIter = config[val_iZone]->GetTimeIter(); - bool StopCalc; - unsigned long iIncrement; - unsigned long nIncrements = config[val_iZone]->GetNumberIncrements(); - bool nonlinear = (config[val_iZone]->GetGeometricConditions() == LARGE_DEFORMATIONS); - bool linear = (config[val_iZone]->GetGeometricConditions() == SMALL_DEFORMATIONS); - bool disc_adj_fem = (config[val_iZone]->GetKind_Solver() == DISC_ADJ_FEM); - bool incremental_load = config[val_iZone]->GetIncrementalLoad(); // Loads applied in steps + const unsigned long TimeIter = config[val_iZone]->GetTimeIter(); + const unsigned long nIncrements = config[val_iZone]->GetNumberIncrements(); + + const bool nonlinear = (config[val_iZone]->GetGeometricConditions() == LARGE_DEFORMATIONS); + const bool linear = (config[val_iZone]->GetGeometricConditions() == SMALL_DEFORMATIONS); + const bool disc_adj_fem = config[val_iZone]->GetDiscrete_Adjoint(); + + /*--- Loads applied in steps (not used for discrete adjoint). ---*/ + const bool incremental_load = config[val_iZone]->GetIncrementalLoad() && !disc_adj_fem; CIntegration* feaIntegration = integration[val_iZone][val_iInst][FEA_SOL]; CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL]; @@ -1234,11 +1235,10 @@ void CFEAIteration::Iterate(COutput *output, config[val_iZone]->SetGlobalParam(FEM_ELASTICITY, RUNTIME_FEA_SYS); if (linear) { + /*--- Run the (one) iteration ---*/ config[val_iZone]->SetInnerIter(0); - /*--- Run the iteration ---*/ - feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, val_iInst); @@ -1285,16 +1285,10 @@ void CFEAIteration::Iterate(COutput *output, /*--- THIS IS THE INCREMENTAL LOAD APPROACH (only makes sense for nonlinear) ---*/ - /*--- Set the initial condition: store the current solution as Solution_Old ---*/ - - feaSolver->SetInitialCondition(geometry[val_iZone][val_iInst], - solver[val_iZone][val_iInst], config[val_iZone], TimeIter); - /*--- Assume the initial load increment as 1.0 ---*/ - loadIncrement = 1.0; - feaSolver->SetLoad_Increment(loadIncrement); - feaSolver->SetForceCoeff(loadIncrement); + feaSolver->SetLoad_Increment(0, 1.0); + feaSolver->SetForceCoeff(1.0); /*--- Run two nonlinear iterations to check if incremental loading can be skipped ---*/ @@ -1302,34 +1296,23 @@ void CFEAIteration::Iterate(COutput *output, config[val_iZone]->SetInnerIter(IntIter); - /*--- Run the first iteration ---*/ - feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, val_iInst); - /*--- Write the convergence history (first computes Von Mises stress) ---*/ - - 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); } - bool meetCriteria; - su2double Residual_UTOL, Residual_RTOL, Residual_ETOL; - su2double Criteria_UTOL, Criteria_RTOL, Criteria_ETOL; + /*--- Early return if we already meet the convergence criteria. ---*/ + if (StopCalc) return; - Criteria_UTOL = config[val_iZone]->GetIncLoad_Criteria(0); - Criteria_RTOL = config[val_iZone]->GetIncLoad_Criteria(1); - Criteria_ETOL = config[val_iZone]->GetIncLoad_Criteria(2); + /*--- Check user-defined criteria to either increment loads or continue with NR iterations. ---*/ - Residual_UTOL = log10(feaSolver->LinSysSol.norm()); - Residual_RTOL = log10(feaSolver->LinSysRes.norm()); - Residual_ETOL = log10(feaSolver->LinSysSol.dot(feaSolver->LinSysRes)); + bool meetCriteria = true; + for (int i = 0; i < 3; ++i) + meetCriteria &= (log10(feaSolver->GetRes_FEM(i)) < config[val_iZone]->GetIncLoad_Criteria(i)); - meetCriteria = ( ( Residual_UTOL < Criteria_UTOL ) && - ( Residual_RTOL < Criteria_RTOL ) && - ( Residual_ETOL < Criteria_ETOL ) ); - - /*--- If the criteria is met, i.e. the load is not "too big", continue the regular calculation ---*/ + /*--- If the criteria is met, i.e. the load is not too large, continue the regular calculation. ---*/ if (meetCriteria) { @@ -1350,30 +1333,31 @@ void CFEAIteration::Iterate(COutput *output, } - /*--- If the criteria is not met, a whole set of subiterations for the different loads must be done ---*/ + /*--- If the criteria is not met, a whole set of subiterations for the different loads must be done. ---*/ else { - /*--- Here we have to restore the solution to the one before testing the criteria ---*/ - /*--- Retrieve the Solution_Old as the current solution before subiterating ---*/ + /*--- Restore solution to initial. Because we ramp the load from zero, in multizone cases it is not + * adequate to take "old values" as those will be for maximum loading on the previous outer iteration. ---*/ - feaSolver->ResetInitialCondition(geometry[val_iZone][val_iInst], - solver[val_iZone][val_iInst], config[val_iZone], TimeIter); + feaSolver->SetInitialCondition(geometry[val_iZone][val_iInst], + solver[val_iZone][val_iInst], + config[val_iZone], TimeIter); /*--- For the number of increments ---*/ - for (iIncrement = 0; iIncrement < nIncrements; iIncrement++) { + for (auto iIncrement = 1ul; iIncrement <= nIncrements; iIncrement++) { /*--- Set the load increment and the initial condition, and output the * parameters of UTOL, RTOL, ETOL for the previous iteration ---*/ - loadIncrement = (iIncrement + 1.0) * (1.0 / nIncrements); - feaSolver->SetLoad_Increment(loadIncrement); + su2double loadIncrement = su2double(iIncrement) / nIncrements; + feaSolver->SetLoad_Increment(iIncrement, loadIncrement); /*--- Set the convergence monitor to false, to force the solver to converge every subiteration ---*/ output->SetConvergence(false); if (rank == MASTER_NODE) - cout << "\nIncremental load: increment " << iIncrement + 1 << endl; + cout << "\nIncremental load: increment " << iIncrement << endl; /*--- Newton-Raphson subiterations ---*/ @@ -1384,7 +1368,6 @@ void CFEAIteration::Iterate(COutput *output, feaIntegration->Structural_Iteration(geometry, solver, numerics, config, RUNTIME_FEA_SYS, val_iZone, val_iInst); - /*--- Write the convergence history (first, compute Von Mises stress) ---*/ StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox, val_iZone, INST_0); @@ -1393,41 +1376,12 @@ void CFEAIteration::Iterate(COutput *output, } } - + /*--- Just to be sure, set default increment settings. ---*/ + feaSolver->SetLoad_Increment(0, 1.0); } } - /*--- Finally, we need to compute the objective function, in case that we are running a discrete adjoint solver... ---*/ - - switch (config[val_iZone]->GetKind_ObjFunc()) { - case REFERENCE_GEOMETRY: - if ((config[val_iZone]->GetDV_FEA() == YOUNG_MODULUS) || (config[val_iZone]->GetDV_FEA() == DENSITY_VAL)) { - feaSolver->Stiffness_Penalty(geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0], - numerics[val_iZone][val_iInst][MESH_0][FEA_SOL], config[val_iZone]); - } - feaSolver->Compute_OFRefGeom(geometry[val_iZone][val_iInst][MESH_0], - solver[val_iZone][val_iInst][MESH_0], config[val_iZone]); - break; - case REFERENCE_NODE: - if ((config[val_iZone]->GetDV_FEA() == YOUNG_MODULUS) || (config[val_iZone]->GetDV_FEA() == DENSITY_VAL)) { - feaSolver->Stiffness_Penalty(geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0], - numerics[val_iZone][val_iInst][MESH_0][FEA_SOL], config[val_iZone]); - } - feaSolver->Compute_OFRefNode(geometry[val_iZone][val_iInst][MESH_0], - solver[val_iZone][val_iInst][MESH_0], config[val_iZone]); - break; - case VOLUME_FRACTION: - case TOPOL_DISCRETENESS: - feaSolver->Compute_OFVolFrac(geometry[val_iZone][val_iInst][MESH_0], - solver[val_iZone][val_iInst][MESH_0], config[val_iZone]); - break; - case TOPOL_COMPLIANCE: - feaSolver->Compute_OFCompliance(geometry[val_iZone][val_iInst][MESH_0], - solver[val_iZone][val_iInst][MESH_0], config[val_iZone]); - break; - } - } void CFEAIteration::Update(COutput *output, @@ -1442,18 +1396,12 @@ void CFEAIteration::Update(COutput *output, unsigned short val_iZone, unsigned short val_iInst) { - su2double Physical_dt, Physical_t; - unsigned long TimeIter = config[val_iZone]->GetTimeIter(); - bool dynamic = (config[val_iZone]->GetTime_Domain()); // Dynamic problems - bool fsi = config[val_iZone]->GetFSI_Simulation(); // Fluid-Structure Interaction problems + const auto TimeIter = config[val_iZone]->GetTimeIter(); + 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]; - /*----------------- Compute averaged nodal stress and reactions ------------------------*/ - - feaSolver->Compute_NodalStress(geometry[val_iZone][val_iInst][MESH_0], - numerics[val_iZone][val_iInst][MESH_0][FEA_SOL], config[val_iZone]); - /*----------------- Update structural solver ----------------------*/ if (dynamic) { @@ -1464,17 +1412,16 @@ void CFEAIteration::Update(COutput *output, /*--- Verify convergence criteria (based on total time) ---*/ - Physical_dt = config[val_iZone]->GetDelta_DynTime(); - Physical_t = (TimeIter+1)*Physical_dt; - if (Physical_t >= config[val_iZone]->GetTotal_DynTime()) + const su2double Physical_dt = config[val_iZone]->GetDelta_DynTime(); + const su2double Physical_t = (TimeIter+1)*Physical_dt; + if (Physical_t >= config[val_iZone]->GetTotal_DynTime()) integration[val_iZone][val_iInst][FEA_SOL]->SetConvergence(true); } 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]->GetKind_TimeIntScheme_FEA() == NEWMARK_IMPLICIT) { - feaSolver->ImplicitNewmark_Relaxation(geometry[val_iZone][val_iInst][MESH_0], - solver[val_iZone][val_iInst][MESH_0], config[val_iZone]); + feaSolver->ImplicitNewmark_Relaxation(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); } } @@ -1494,16 +1441,10 @@ void CFEAIteration::Predictor(COutput *output, CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL]; - /*--- Predict displacements ---*/ - - feaSolver->PredictStruct_Displacement(geometry[val_iZone][val_iInst], config[val_iZone], solver[val_iZone][val_iInst]); - - /*--- For parallel simulations we need to communicate the predicted solution before updating the fluid mesh ---*/ - - feaSolver->InitiateComms(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], SOLUTION_PRED); - feaSolver->CompleteComms(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], SOLUTION_PRED); + feaSolver->PredictStruct_Displacement(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); } + void CFEAIteration::Relaxation(COutput *output, CIntegration ****integration, CGeometry ****geometry, @@ -1522,17 +1463,12 @@ void CFEAIteration::Relaxation(COutput *output, /*------------------- Compute the coefficient ---------------------*/ - feaSolver->ComputeAitken_Coefficient(geometry[val_iZone][INST_0], config[val_iZone], - solver[val_iZone][INST_0], config[val_iZone]->GetOuterIter()); + feaSolver->ComputeAitken_Coefficient(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone], + config[val_iZone]->GetOuterIter()); /*----------------- Set the relaxation parameter ------------------*/ - feaSolver->SetAitken_Relaxation(geometry[val_iZone][INST_0], config[val_iZone], solver[val_iZone][INST_0]); - - /*----------------- Communicate the predicted solution and the old one ------------------*/ - - feaSolver->InitiateComms(geometry[val_iZone][INST_0][MESH_0], config[val_iZone], SOLUTION_PRED_OLD); - feaSolver->CompleteComms(geometry[val_iZone][INST_0][MESH_0], config[val_iZone], SOLUTION_PRED_OLD); + feaSolver->SetAitken_Relaxation(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); } @@ -1555,12 +1491,9 @@ bool CFEAIteration::Monitor(COutput *output, #endif UsedTime = StopTime - StartTime; - solver[val_iZone][val_iInst][MESH_0][FEA_SOL]->Compute_NodalStress(geometry[val_iZone][val_iInst][MESH_0], - numerics[val_iZone][val_iInst][MESH_0][FEA_SOL], - config[val_iZone]); - if (config[val_iZone]->GetMultizone_Problem() || config[val_iZone]->GetSinglezone_Driver()){ - output->SetHistory_Output(geometry[val_iZone][INST_0][MESH_0], solver[val_iZone][INST_0][MESH_0], + output->SetHistory_Output(geometry[val_iZone][val_iInst][MESH_0], + solver[val_iZone][val_iInst][MESH_0], config[val_iZone], config[val_iZone]->GetTimeIter(), config[val_iZone]->GetOuterIter(), config[val_iZone]->GetInnerIter()); } @@ -1596,13 +1529,13 @@ void CFEAIteration::Solve(COutput *output, /*------------------ Structural subiteration ----------------------*/ Iterate(output, integration, geometry, solver, numerics, config, - surface_movement, grid_movement, FFDBox, val_iZone, INST_0); + surface_movement, grid_movement, FFDBox, val_iZone, val_iInst); /*--- Write the convergence history for the structure (only screen output) ---*/ // if (multizone) output->SetConvHistory_Body(geometry, solver, config, integration, false, 0.0, val_iZone, INST_0); /*--- Set the structural convergence to false (to make sure outer subiterations converge) ---*/ - integration[val_iZone][INST_0][FEA_SOL]->SetConvergence(false); + integration[val_iZone][val_iInst][FEA_SOL]->SetConvergence(false); } @@ -2115,57 +2048,44 @@ void CDiscAdjFluidIteration::Iterate(COutput *output, CSurfaceMovement **surface_movement, CVolumetricMovement ***volume_grid_movement, CFreeFormDefBox*** FFDBox, - unsigned short val_iZone, - unsigned short val_iInst) { - - unsigned short Kind_Solver = config[val_iZone]->GetKind_Solver(); - bool frozen_visc = config[val_iZone]->GetFrozen_Visc_Disc(); - bool heat = config[val_iZone]->GetWeakly_Coupled_Heat(); + unsigned short iZone, + unsigned short iInst) { + bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); + bool heat = config[iZone]->GetWeakly_Coupled_Heat(); /*--- Extract the adjoints of the conservative input variables and store them for the next iteration ---*/ - if ((Kind_Solver == DISC_ADJ_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_RANS) || (Kind_Solver == DISC_ADJ_EULER) || - (Kind_Solver == DISC_ADJ_INC_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_INC_RANS) || (Kind_Solver == DISC_ADJ_INC_EULER)) { - - solver[val_iZone][val_iInst][MESH_0][ADJFLOW_SOL]->ExtractAdjoint_Solution(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); - - solver[val_iZone][val_iInst][MESH_0][ADJFLOW_SOL]->ExtractAdjoint_Variables(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); + if (config[iZone]->GetFluidProblem()) { + solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone]); + solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]); } if (turbulent && !frozen_visc) { - solver[val_iZone][val_iInst][MESH_0][ADJTURB_SOL]->ExtractAdjoint_Solution(geometry[val_iZone][val_iInst][MESH_0], - config[val_iZone]); + solver[iZone][iInst][MESH_0][ADJTURB_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone]); } if (heat) { - - solver[val_iZone][val_iInst][MESH_0][ADJHEAT_SOL]->ExtractAdjoint_Solution(geometry[val_iZone][val_iInst][MESH_0], - config[val_iZone]); + solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone]); } - if (config[val_iZone]->AddRadiation()) { - - solver[val_iZone][val_iInst][MESH_0][ADJRAD_SOL]->ExtractAdjoint_Solution(geometry[val_iZone][val_iInst][MESH_0], - config[val_iZone]); - - solver[val_iZone][val_iInst][MESH_0][ADJRAD_SOL]->ExtractAdjoint_Variables(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); + if (config[iZone]->AddRadiation()) { + solver[iZone][iInst][MESH_0][ADJRAD_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone]); + solver[iZone][iInst][MESH_0][ADJRAD_SOL]->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]); } + } void CDiscAdjFluidIteration::InitializeAdjoint(CSolver *****solver, CGeometry ****geometry, CConfig **config, unsigned short iZone, unsigned short iInst){ - unsigned short Kind_Solver = config[iZone]->GetKind_Solver(); bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); bool heat = config[iZone]->GetWeakly_Coupled_Heat(); bool interface_boundary = (config[iZone]->GetnMarker_Fluid_Load() > 0); /*--- Initialize the adjoints the conservative variables ---*/ - if ((Kind_Solver == DISC_ADJ_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_RANS) || (Kind_Solver == DISC_ADJ_EULER) || - (Kind_Solver == DISC_ADJ_INC_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_INC_RANS) || (Kind_Solver == DISC_ADJ_INC_EULER)) { - + if (config[iZone]->GetFluidProblem()) { solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]); } @@ -2190,17 +2110,14 @@ void CDiscAdjFluidIteration::InitializeAdjoint(CSolver *****solver, CGeometry ** void CDiscAdjFluidIteration::RegisterInput(CSolver *****solver, CGeometry ****geometry, CConfig **config, unsigned short iZone, unsigned short iInst, unsigned short kind_recording){ - unsigned short Kind_Solver = config[iZone]->GetKind_Solver(); bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); bool heat = config[iZone]->GetWeakly_Coupled_Heat(); - if (kind_recording == FLOW_CONS_VARS || kind_recording == COMBINED){ + if (kind_recording == SOLUTION_VARIABLES || kind_recording == SOLUTION_AND_MESH) { /*--- Register flow and turbulent variables as input ---*/ - if ((Kind_Solver == DISC_ADJ_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_RANS) || (Kind_Solver == DISC_ADJ_EULER) || - (Kind_Solver == DISC_ADJ_INC_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_INC_RANS) || (Kind_Solver == DISC_ADJ_INC_EULER)) { - + if (config[iZone]->GetFluidProblem()) { solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]); solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]); @@ -2218,29 +2135,8 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver *****solver, CGeometry ****ge solver[iZone][iInst][MESH_0][ADJRAD_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]); } } - if (kind_recording == MESH_COORDS){ - - /*--- Register node coordinates as input ---*/ - - geometry[iZone][iInst][MESH_0]->RegisterCoordinates(config[iZone]); - - } - - if (kind_recording == FLOW_CROSS_TERM){ - - /*--- Register flow and turbulent variables as input ---*/ - - solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]); - - if (turbulent && !frozen_visc){ - solver[iZone][iInst][MESH_0][ADJTURB_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]); - } - if (config[iZone]->AddRadiation()) { - solver[iZone][iInst][MESH_0][ADJRAD_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]); - } - } - if (kind_recording == GEOMETRY_CROSS_TERM){ + if (kind_recording == MESH_COORDS){ /*--- Register node coordinates as input ---*/ @@ -2264,27 +2160,29 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver *****solver, CGeometry ****ge void CDiscAdjFluidIteration::SetRecording(CSolver *****solver, CGeometry ****geometry, CConfig **config, - unsigned short val_iZone, - unsigned short val_iInst, + unsigned short iZone, + unsigned short iInst, unsigned short kind_recording) { - unsigned short iMesh; + bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); /*--- Prepare for recording by resetting the solution to the initial converged solution ---*/ - solver[val_iZone][val_iInst][MESH_0][ADJFEA_SOL]->SetRecording(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); + if (solver[iZone][iInst][MESH_0][ADJFEA_SOL]) { + solver[iZone][iInst][MESH_0][ADJFEA_SOL]->SetRecording(geometry[iZone][iInst][MESH_0], config[iZone]); + } - for (iMesh = 0; iMesh <= config[val_iZone]->GetnMGLevels(); iMesh++){ - solver[val_iZone][val_iInst][iMesh][ADJFLOW_SOL]->SetRecording(geometry[val_iZone][val_iInst][iMesh], config[val_iZone]); + for (auto iMesh = 0u; iMesh <= config[iZone]->GetnMGLevels(); iMesh++){ + solver[iZone][iInst][iMesh][ADJFLOW_SOL]->SetRecording(geometry[iZone][iInst][iMesh], config[iZone]); } - if ((config[val_iZone]->GetKind_Solver() == DISC_ADJ_RANS || config[val_iZone]->GetKind_Solver() == DISC_ADJ_INC_RANS) && !config[val_iZone]->GetFrozen_Visc_Disc()) { - solver[val_iZone][val_iInst][MESH_0][ADJTURB_SOL]->SetRecording(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); + if (turbulent && !frozen_visc) { + solver[iZone][iInst][MESH_0][ADJTURB_SOL]->SetRecording(geometry[iZone][iInst][MESH_0], config[iZone]); } - if (config[val_iZone]->GetWeakly_Coupled_Heat()) { - solver[val_iZone][val_iInst][MESH_0][ADJHEAT_SOL]->SetRecording(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); + if (config[iZone]->GetWeakly_Coupled_Heat()) { + solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->SetRecording(geometry[iZone][iInst][MESH_0], config[iZone]); } - if (config[val_iZone]->GetKind_RadiationModel() != NONE) { - solver[val_iZone][INST_0][MESH_0][ADJRAD_SOL]->SetRecording(geometry[val_iZone][INST_0][MESH_0], config[val_iZone]); + if (config[iZone]->AddRadiation()) { + solver[iZone][INST_0][MESH_0][ADJRAD_SOL]->SetRecording(geometry[iZone][INST_0][MESH_0], config[iZone]); } } @@ -2299,8 +2197,8 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver *****solver, bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); bool heat = config[iZone]->GetWeakly_Coupled_Heat(); - if ((kind_recording == MESH_COORDS) || (kind_recording == NONE) || - (kind_recording == GEOMETRY_CROSS_TERM) || (kind_recording == ALL_VARIABLES)){ + + if ((kind_recording == MESH_COORDS) || (kind_recording == NONE) || (kind_recording == SOLUTION_AND_MESH)){ /*--- Update geometry to get the influence on other geometry variables (normals, volume etc) ---*/ @@ -2337,26 +2235,20 @@ void CDiscAdjFluidIteration::SetDependencies(CSolver *****solver, void CDiscAdjFluidIteration::RegisterOutput(CSolver *****solver, CGeometry ****geometry, CConfig **config, COutput* output, unsigned short iZone, unsigned short iInst){ - unsigned short Kind_Solver = config[iZone]->GetKind_Solver(); bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); bool heat = config[iZone]->GetWeakly_Coupled_Heat(); bool interface_boundary = (config[iZone]->GetnMarker_Fluid_Load() > 0); - if ((Kind_Solver == DISC_ADJ_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_RANS) || (Kind_Solver == DISC_ADJ_EULER) || - (Kind_Solver == DISC_ADJ_INC_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_INC_RANS) || (Kind_Solver == DISC_ADJ_INC_EULER)) { - /*--- Register conservative variables as output of the iteration ---*/ - solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0],config[iZone]); - + if (config[iZone]->GetFluidProblem()){ + solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); } if (turbulent && !frozen_visc){ - solver[iZone][iInst][MESH_0][ADJTURB_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], - config[iZone]); + solver[iZone][iInst][MESH_0][ADJTURB_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); } if (heat){ - solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], - config[iZone]); + solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); } if (config[iZone]->AddRadiation()){ solver[iZone][iInst][MESH_0][ADJRAD_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); @@ -2366,35 +2258,6 @@ void CDiscAdjFluidIteration::RegisterOutput(CSolver *****solver, CGeometry ****g } } -void CDiscAdjFluidIteration::InitializeAdjoint_CrossTerm(CSolver *****solver, CGeometry ****geometry, CConfig **config, unsigned short iZone, unsigned short iInst){ - - unsigned short Kind_Solver = config[iZone]->GetKind_Solver(); - bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); - - /*--- Initialize the adjoint of the objective function (typically with 1.0) ---*/ - - solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->SetAdj_ObjFunc(geometry[iZone][iInst][MESH_0], config[iZone]); - - /*--- Initialize the adjoints the conservative variables ---*/ - - if ((Kind_Solver == DISC_ADJ_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_RANS) || (Kind_Solver == DISC_ADJ_EULER) || - (Kind_Solver == DISC_ADJ_INC_NAVIER_STOKES) || (Kind_Solver == DISC_ADJ_INC_RANS) || (Kind_Solver == DISC_ADJ_INC_EULER)) { - - solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], - config[iZone]); -} - - if (turbulent && !frozen_visc) { - solver[iZone][iInst][MESH_0][ADJTURB_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], - config[iZone]); - } - - if (config[iZone]->AddRadiation()) { - solver[iZone][iInst][MESH_0][ADJRAD_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]); - } -} - - void CDiscAdjFluidIteration::Update(COutput *output, CIntegration ****integration, CGeometry ****geometry, @@ -2419,6 +2282,7 @@ void CDiscAdjFluidIteration::Update(COutput *output, } } } + bool CDiscAdjFluidIteration::Monitor(COutput *output, CIntegration ****integration, CGeometry ****geometry, @@ -2670,7 +2534,7 @@ void CDiscAdjFEAIteration::SetRecording(COutput *output, /*--- Clear indices of coupling variables ---*/ - SetDependencies(solver, geometry, numerics, config, val_iZone, val_iInst, ALL_VARIABLES); + SetDependencies(solver, geometry, numerics, config, val_iZone, val_iInst, SOLUTION_AND_MESH); /*--- Run one iteration while tape is passive - this clears all indices ---*/ @@ -2751,9 +2615,6 @@ void CDiscAdjFEAIteration::RegisterInput(CSolver *****solver, CGeometry ****geom /*--- Register variables as input ---*/ solver[iZone][iInst][MESH_0][ADJFEA_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]); - - /*--- Both need to be registered regardless of kind_recording for structural shape derivatives to work properly. - Otherwise, the code simply diverges as the FEM_CROSS_TERM_GEOMETRY breaks! (no idea why) for this term we register but do not extract! ---*/ } else { /*--- Register topology optimization densities (note direct solver) ---*/ @@ -2850,7 +2711,7 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver *****solver, CGeometry ****ge /*--- FSI specific dependencies. ---*/ if (fsi) { /*--- Set relation between solution and predicted displacements, which are the transferred ones. ---*/ - dir_solver->PredictStruct_Displacement(nullptr, config[iZone], solver[iZone][iInst]); + dir_solver->PredictStruct_Displacement(structural_geometry, config[iZone]); } /*--- MPI dependencies. ---*/ @@ -2896,19 +2757,6 @@ void CDiscAdjFEAIteration::InitializeAdjoint(CSolver *****solver, CGeometry **** } - -void CDiscAdjFEAIteration::InitializeAdjoint_CrossTerm(CSolver *****solver, CGeometry ****geometry, CConfig **config, unsigned short iZone, unsigned short iInst){ - - /*--- Initialize the adjoint of the objective function (typically with 1.0) ---*/ - - solver[iZone][iInst][MESH_0][ADJFEA_SOL]->SetAdj_ObjFunc(geometry[iZone][iInst][MESH_0], config[iZone]); - - /*--- Initialize the adjoints the conservative variables ---*/ - - solver[iZone][iInst][MESH_0][ADJFEA_SOL]->SetAdjoint_Output(geometry[iZone][iInst][MESH_0], config[iZone]); - -} - void CDiscAdjFEAIteration::Update(COutput *output, CIntegration ****integration, CGeometry ****geometry, @@ -3270,7 +3118,7 @@ void CDiscAdjHeatIteration::RegisterInput(CSolver *****solver, unsigned short iZone, unsigned short iInst, unsigned short kind_recording){ - if (kind_recording == FLOW_CONS_VARS || kind_recording == COMBINED){ + if (kind_recording == SOLUTION_VARIABLES || kind_recording == SOLUTION_AND_MESH){ /*--- Register flow and turbulent variables as input ---*/ @@ -3295,8 +3143,7 @@ void CDiscAdjHeatIteration::SetDependencies(CSolver *****solver, unsigned short iZone, unsigned short iInst, unsigned short kind_recording){ - if ((kind_recording == MESH_COORDS) || (kind_recording == NONE) || - (kind_recording == GEOMETRY_CROSS_TERM) || (kind_recording == ALL_VARIABLES)){ + if ((kind_recording == MESH_COORDS) || (kind_recording == NONE) || (kind_recording == SOLUTION_AND_MESH)){ /*--- Update geometry to get the influence on other geometry variables (normals, volume etc) ---*/ diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp index 64109badd937..8999609493a8 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp @@ -187,9 +187,12 @@ void CFEAElasticity::Compute_Mass_Matrix(CElement *element, const CConfig *confi unsigned short iGauss, nGauss; unsigned short iNode, jNode, nNode; - su2double Weight, Jac_X; + su2double Weight, Jac_X, val_Mab; - su2double val_Mab; + /*--- Register pre-accumulation inputs, density and reference coords. ---*/ + AD::StartPreacc(); + AD::SetPreaccIn(Rho_s); + element->SetPreaccIn_Coords(false); element->ClearElement(); /*--- Restarts the element: avoids adding over previous results in other elements --*/ element->ComputeGrad_Linear(); /*--- Need to compute the gradients to obtain the Jacobian ---*/ @@ -227,6 +230,10 @@ void CFEAElasticity::Compute_Mass_Matrix(CElement *element, const CConfig *confi } + /*--- Register the mass matrix as preaccumulation output. ---*/ + element->SetPreaccOut_Mab(); + AD::EndPreacc(); + } diff --git a/SU2_CFD/src/output/CAdjElasticityOutput.cpp b/SU2_CFD/src/output/CAdjElasticityOutput.cpp index ac50dfcb397b..01ac0a3ecb3b 100644 --- a/SU2_CFD/src/output/CAdjElasticityOutput.cpp +++ b/SU2_CFD/src/output/CAdjElasticityOutput.cpp @@ -101,6 +101,8 @@ void CAdjElasticityOutput::SetHistoryOutputFields(CConfig *config){ AddHistoryOutput("SENS_E", "Sens[E]", ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", ""); AddHistoryOutput("SENS_NU","Sens[Nu]", ScreenOutputFormat::SCIENTIFIC, "SENSITIVITY", ""); + AddHistoryOutput("COMBO", "ObjFun", ScreenOutputFormat::SCIENTIFIC, "COMBO", "", HistoryFieldType::COEFFICIENT); + } inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) { @@ -110,23 +112,25 @@ inline void CAdjElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *ge 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){ + if (config->GetnElasticityMod() == 1) { Total_SensE = solver[ADJFEA_SOL]->GetGlobal_Sens_E(0); Total_SensNu = solver[ADJFEA_SOL]->GetGlobal_Sens_Nu(0); } - else{ - // TODO: Update this and change tests + else { for (unsigned short iVar = 0; iVar < config->GetnElasticityMod(); iVar++){ - Total_SensE += pow(solver[ADJFEA_SOL]->GetGlobal_Sens_E(0),2); - Total_SensNu += pow(solver[ADJFEA_SOL]->GetGlobal_Sens_Nu(0),2); + Total_SensE += pow(solver[ADJFEA_SOL]->GetGlobal_Sens_E(iVar),2); + Total_SensNu += pow(solver[ADJFEA_SOL]->GetGlobal_Sens_Nu(iVar),2); } - Total_SensE = sqrt(Total_SensE); - Total_SensNu = sqrt(Total_SensNu); + Total_SensE = sqrt(Total_SensE); + Total_SensNu = sqrt(Total_SensNu); } SetHistoryOutputValue("SENS_E", Total_SensE); SetHistoryOutputValue("SENS_NU", Total_SensNu); + SetHistoryOutputValue("COMBO", solver[FEA_SOL]->GetTotal_ComboObj()); + } void CAdjElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){ diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index 0394a6e1eb66..41521088d3c2 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -115,10 +115,9 @@ void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS SetHistoryOutputValue("RMS_DISP_Z", log10(fea_solver->GetRes_RMS(2))); } } else if (nonlinear_analysis){ - SetHistoryOutputValue("RMS_UTOL", log10(fea_solver->LinSysSol.norm())); - SetHistoryOutputValue("RMS_RTOL", log10(fea_solver->LinSysRes.norm())); - SetHistoryOutputValue("RMS_ETOL", log10(fea_solver->LinSysSol.dot(fea_solver->LinSysRes))); - + SetHistoryOutputValue("RMS_UTOL", log10(fea_solver->GetRes_FEM(0))); + SetHistoryOutputValue("RMS_RTOL", log10(fea_solver->GetRes_FEM(1))); + SetHistoryOutputValue("RMS_ETOL", log10(fea_solver->GetRes_FEM(2))); } if (multiZone){ @@ -134,6 +133,8 @@ void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS SetHistoryOutputValue("LINSOL_ITER", fea_solver->GetIterLinSolver()); SetHistoryOutputValue("LINSOL_RESIDUAL", log10(fea_solver->GetResLinSolver())); + SetHistoryOutputValue("COMBO", fea_solver->GetTotal_ComboObj()); + } void CElasticityOutput::SetHistoryOutputFields(CConfig *config){ @@ -159,6 +160,8 @@ void CElasticityOutput::SetHistoryOutputFields(CConfig *config){ AddHistoryOutput("LOAD_INCREMENT", "Load[%]", ScreenOutputFormat::PERCENT, "", "LOAD_INCREMENT"); AddHistoryOutput("LOAD_RAMP", "Load_Ramp", ScreenOutputFormat::FIXED, "", "LOAD_RAMP"); + AddHistoryOutput("COMBO", "ObjFun", ScreenOutputFormat::SCIENTIFIC, "COMBO", "", HistoryFieldType::COEFFICIENT); + } void CElasticityOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolver **solver, unsigned long iPoint){ @@ -237,5 +240,3 @@ bool CElasticityOutput::SetInit_Residuals(CConfig *config){ return (config->GetTime_Domain() == NO && (curInnerIter == 0)); } - - diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 5a24ace56ff1..1c1864f5f667 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -48,6 +48,7 @@ CFEASolver::CFEASolver(bool mesh_deform_mode) : CSolver(mesh_deform_mode) { Total_CFEA = 0.0; WAitken_Dyn = 0.0; WAitken_Dyn_tn1 = 0.0; + idxIncrement = 0; loadIncrement = 1.0; element_container = new CElement** [MAX_TERMS](); @@ -123,6 +124,7 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CSolver() { Total_CFEA = 0.0; WAitken_Dyn = 0.0; WAitken_Dyn_tn1 = 0.0; + idxIncrement = 0; loadIncrement = 0.0; SetFSI_ConvValue(0,0.0); @@ -188,9 +190,6 @@ CFEASolver::CFEASolver(CGeometry *geometry, CConfig *config) : CSolver() { /*--- Initialization of linear solver structures ---*/ LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); - - LinSysAux.Initialize(nPoint, nPointDomain, nVar, 0.0); - LinSysReact.Initialize(nPoint, nPointDomain, nVar, 0.0); /*--- Initialize structures for hybrid-parallel mode. ---*/ @@ -530,7 +529,7 @@ void CFEASolver::Set_Prestretch(CGeometry *geometry, CConfig *config) { iPoint_Local = Global2Local[iPoint_Global]; - su2double Sol[3] = {0.0, 0.0, 0.0}; + su2double Sol[MAXNVAR] = {0.0}; if (nDim == 2) point_line >> Sol[0] >> Sol[1] >> index; if (nDim == 3) point_line >> Sol[0] >> Sol[1] >> Sol[2] >> index; @@ -683,7 +682,7 @@ void CFEASolver::Set_ReferenceGeometry(CGeometry *geometry, CConfig *config) { if (iPoint_Local >= 0) { - su2double Sol[3] = {0.0, 0.0, 0.0}; + su2double Sol[MAXNVAR] = {0.0}; if (nDim == 2){ Sol[0] = PrintingToolbox::stod(point_line[3]); @@ -826,7 +825,7 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, /* * FSI loads (computed upstream) need to be integrated if a nonconservative interpolation scheme is in use */ - if (fsi && first_iter) Integrate_FSI_Loads(geometry, config); + if (fsi && first_iter && (idxIncrement==0)) Integrate_FSI_Loads(geometry, config); /*--- Next call to Preprocessing will not be "initial_calc" and linear operations will not be repeated. ---*/ initial_calc = false; @@ -835,21 +834,23 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, void CFEASolver::SetInitialCondition(CGeometry **geometry, CSolver ***solver_container, CConfig *config, unsigned long TimeIter) { - /*--- We store the current solution as "Solution Old", for the case that we need to retrieve it ---*/ - - if (config->GetIncrementalLoad()) nodes->Set_OldSolution(); - -} - -void CFEASolver::ResetInitialCondition(CGeometry **geometry, CSolver ***solver_container, CConfig *config, unsigned long TimeIter) { - - /*--- We store the current solution as "Solution Old", for the case that we need to retrieve it ---*/ - - if (config->GetIncrementalLoad()) nodes->Set_Solution(); - + SU2_OMP_PARALLEL + { + 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); + } + else { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) + nodes->SetSolution(iPoint, nodes->GetPrestretch(iPoint)); + } + } // end parallel } -void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, CConfig *config) { +void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { const bool topology_mode = config->GetTopology_Optimization(); const su2double simp_exponent = config->GetSIMP_Exponent(); @@ -884,7 +885,7 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, CElement* element = element_container[FEA_TERM][EL_KIND+thread*MAX_FE_KINDS]; /*--- For the number of nodes, get the coordinates and cache the point indices. ---*/ - unsigned long indexNode[MAXNNODE]; + unsigned long indexNode[MAXNNODE_3D]; for (iNode = 0; iNode < nNodes; iNode++) { @@ -938,7 +939,7 @@ void CFEASolver::Compute_StiffMatrix(CGeometry *geometry, CNumerics **numerics, } -void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumerics **numerics, CConfig *config) { +void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { const bool prestretch_fem = config->GetPrestretch(); const bool de_effects = config->GetDE_Effects(); @@ -977,7 +978,7 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri CElement* de_elem = element_container[DE_TERM][EL_KIND+thread*MAX_FE_KINDS]; /*--- For the number of nodes, we get the coordinates from the connectivity matrix ---*/ - unsigned long indexNode[MAXNNODE]; + unsigned long indexNode[MAXNNODE_3D]; for (iNode = 0; iNode < nNodes; iNode++) { @@ -1078,11 +1079,15 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri } -void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, CConfig *config) { +void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { const bool topology_mode = config->GetTopology_Optimization(); const su2double simp_minstiff = config->GetSIMP_MinStiffness(); + /*--- Never record this method as the mass matrix is passive (but the mass residual is not). ---*/ + const bool ActiveTape = AD::TapeActive(); + AD::StopRecording(); + /*--- Start OpenMP parallel region. ---*/ SU2_OMP_PARALLEL @@ -1111,7 +1116,7 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, C CElement* element = element_container[FEA_TERM][EL_KIND+thread*MAX_FE_KINDS]; /*--- For the number of nodes, get the coordinates and cache the point indices. ---*/ - unsigned long indexNode[MAXNNODE]; + unsigned long indexNode[MAXNNODE_3D]; for (iNode = 0; iNode < nNodes; iNode++) { indexNode[iNode] = geometry->elem[iElem]->GetNode(iNode); @@ -1141,10 +1146,10 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, C for (jNode = 0; jNode < nNodes; jNode++) { auto Mij = MassMatrix.GetBlock(indexNode[iNode], indexNode[jNode]); - su2double Mab = element->Get_Mab(iNode, jNode); + su2double Mab = simp_penalty * element->Get_Mab(iNode, jNode); for (iVar = 0; iVar < nVar; iVar++) - Mij[iVar*(nVar+1)] += simp_penalty*Mab; + Mij[iVar*(nVar+1)] += SU2_TYPE::GetValue(Mab); } if (LockStrategy) omp_unset_lock(&UpdateLocks[indexNode[iNode]]); @@ -1156,90 +1161,86 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, C } // end SU2_OMP_PARALLEL + if (ActiveTape) AD::StartRecording(); + } -void CFEASolver::Compute_MassRes(CGeometry *geometry, CNumerics **numerics, CConfig *config) { +void CFEASolver::Compute_MassRes(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { const bool topology_mode = config->GetTopology_Optimization(); const su2double simp_minstiff = config->GetSIMP_MinStiffness(); - /*--- Start OpenMP parallel region. ---*/ - - SU2_OMP_PARALLEL - { - /*--- Clear vector before calculation. ---*/ - TimeRes.SetValZero(); - SU2_OMP_BARRIER + /*--- Clear vector before calculation. ---*/ + TimeRes.SetValZero(); + SU2_OMP_BARRIER - for(auto color : ElemColoring) { + for(auto color : ElemColoring) { - /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ - SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) - for(auto k = 0ul; k < color.size; ++k) { + /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/ + SU2_OMP_FOR_DYN(nextMultiple(OMP_MIN_SIZE, color.groupSize)) + for(auto k = 0ul; k < color.size; ++k) { - auto iElem = color.indices[k]; + auto iElem = color.indices[k]; - unsigned short iNode, jNode, iDim, iVar; + unsigned short iNode, jNode, iDim, iVar; - int thread = omp_get_thread_num(); + int thread = omp_get_thread_num(); - /*--- Convert VTK type to index in the element container. ---*/ - int EL_KIND; - unsigned short nNodes; - GetElemKindAndNumNodes(geometry->elem[iElem]->GetVTK_Type(), EL_KIND, nNodes); + /*--- Convert VTK type to index in the element container. ---*/ + int EL_KIND; + unsigned short nNodes; + GetElemKindAndNumNodes(geometry->elem[iElem]->GetVTK_Type(), EL_KIND, nNodes); - /*--- Each thread needs a dedicated element. ---*/ - CElement* element = element_container[FEA_TERM][EL_KIND+thread*MAX_FE_KINDS]; + /*--- Each thread needs a dedicated element. ---*/ + CElement* element = element_container[FEA_TERM][EL_KIND+thread*MAX_FE_KINDS]; - /*--- For the number of nodes, get the coordinates and cache the point indices. ---*/ - unsigned long indexNode[MAXNNODE]; + /*--- For the number of nodes, get the coordinates and cache the point indices. ---*/ + unsigned long indexNode[MAXNNODE_3D]; - for (iNode = 0; iNode < nNodes; iNode++) { - indexNode[iNode] = geometry->elem[iElem]->GetNode(iNode); - for (iDim = 0; iDim < nDim; iDim++) { - su2double val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim); - element->SetRef_Coord(iNode, iDim, val_Coord); - } + for (iNode = 0; iNode < nNodes; iNode++) { + indexNode[iNode] = geometry->elem[iElem]->GetNode(iNode); + for (iDim = 0; iDim < nDim; iDim++) { + su2double val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim); + element->SetRef_Coord(iNode, iDim, val_Coord); } + } - /*--- In topology mode determine the penalty to apply to the mass, - * linear function of the physical density. ---*/ - su2double simp_penalty = 1.0; - if (topology_mode) { - simp_penalty = simp_minstiff+(1.0-simp_minstiff)*element_properties[iElem]->GetPhysicalDensity(); - } + /*--- In topology mode determine the penalty to apply to the mass, + * linear function of the physical density. ---*/ + su2double simp_penalty = 1.0; + if (topology_mode) { + simp_penalty = simp_minstiff+(1.0-simp_minstiff)*element_properties[iElem]->GetPhysicalDensity(); + } - /*--- Set the properties of the element and compute its mass matrix. ---*/ - element->Set_ElProperties(element_properties[iElem]); + /*--- Set the properties of the element and compute its mass matrix. ---*/ + element->Set_ElProperties(element_properties[iElem]); - numerics[FEA_TERM + thread*MAX_TERMS]->Compute_Mass_Matrix(element, config); + numerics[FEA_TERM + thread*MAX_TERMS]->Compute_Mass_Matrix(element, config); - /*--- Add contributions of this element to the mass matrix. - * In case we need to use locks we guard the entire update. ---*/ - if (LockStrategy) omp_set_lock(&UpdateLocks[0]); + /*--- Add contributions of this element to the mass residual. + * Equiv. to a matrix vector product with the mass matrix. ---*/ + for (iNode = 0; iNode < nNodes; iNode++) { - for (iNode = 0; iNode < nNodes; iNode++) { - for (jNode = 0; jNode < nNodes; jNode++) { + if (LockStrategy) omp_set_lock(&UpdateLocks[indexNode[iNode]]); - su2double Mab = simp_penalty * element->Get_Mab(iNode, jNode); + for (jNode = 0; jNode < nNodes; jNode++) { - for (iVar = 0; iVar < nVar; iVar++) { - TimeRes[indexNode[iNode]*nVar+iVar] += Mab * TimeRes_Aux(indexNode[iNode],iVar); - TimeRes[indexNode[jNode]*nVar+iVar] += Mab * TimeRes_Aux(indexNode[jNode],iVar); - } - } + su2double Mab = simp_penalty * element->Get_Mab(iNode, jNode); + + for (iVar = 0; iVar < nVar; iVar++) + TimeRes(indexNode[iNode], iVar) += Mab * TimeRes_Aux(indexNode[jNode], iVar); } - if (LockStrategy) omp_unset_lock(&UpdateLocks[0]); - } // end iElem loop + if (LockStrategy) omp_unset_lock(&UpdateLocks[indexNode[iNode]]); + } - } // end color loop + } // end iElem loop - } // end SU2_OMP_PARALLEL + } // end color loop } -void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numerics, CConfig *config) { +void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { const bool prestretch_fem = config->GetPrestretch(); @@ -1276,7 +1277,7 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric CElement* element = element_container[FEA_TERM][EL_KIND+thread*MAX_FE_KINDS]; /*--- For the number of nodes, we get the coordinates from the connectivity matrix ---*/ - unsigned long indexNode[MAXNNODE]; + unsigned long indexNode[MAXNNODE_3D]; for (iNode = 0; iNode < nNodes; iNode++) { @@ -1330,7 +1331,11 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric } -void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, CConfig *config) { +void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { + + /*--- Never record this method as atm it is not differentiable. ---*/ + const bool ActiveTape = AD::TapeActive(); + AD::StopRecording(); const bool prestretch_fem = config->GetPrestretch(); @@ -1377,7 +1382,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, CElement* element = element_container[FEA_TERM][EL_KIND+thread*MAX_FE_KINDS]; /*--- For the number of nodes, we get the coordinates from the connectivity matrix ---*/ - unsigned long indexNode[MAXNNODE]; + unsigned long indexNode[MAXNNODE_3D]; for (iNode = 0; iNode < nNodes; iNode++) { @@ -1558,7 +1563,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, } /*--- Once computed, compute M*TimeRes_Aux ---*/ - MassMatrix.MatrixVectorProduct(TimeRes_Aux,TimeRes,geometry,config); + Compute_MassRes(geometry, numerics, config); /*--- Loop over all the markers ---*/ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) @@ -1609,9 +1614,11 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, } + if (ActiveTape) AD::StartRecording(); + } -void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, CConfig *config) { +void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { /*--- Start OpenMP parallel region. ---*/ @@ -1643,7 +1650,7 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, CCo CElement* element = element_container[FEA_TERM][EL_KIND+thread*MAX_FE_KINDS]; /*--- For the number of nodes, get the coordinates and cache the point indices. ---*/ - unsigned long indexNode[MAXNNODE]; + unsigned long indexNode[MAXNNODE_3D]; for (iNode = 0; iNode < nNodes; iNode++) { indexNode[iNode] = geometry->elem[iElem]->GetNode(iNode); @@ -1686,7 +1693,7 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, CCo } -void CFEASolver::Compute_IntegrationConstants(CConfig *config) { +void CFEASolver::Compute_IntegrationConstants(const CConfig *config) { su2double Delta_t= config->GetDelta_DynTime(); @@ -1737,10 +1744,10 @@ void CFEASolver::Compute_IntegrationConstants(CConfig *config) { } -void CFEASolver::BC_Clamped(CGeometry *geometry, CNumerics *numerics, CConfig *config, unsigned short val_marker) { +void CFEASolver::BC_Clamped(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker) { const bool dynamic = config->GetTime_Domain(); - const su2double zeros[3] = {0.0, 0.0, 0.0}; + const su2double zeros[MAXNVAR] = {0.0}; for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { @@ -1769,11 +1776,11 @@ void CFEASolver::BC_Clamped(CGeometry *geometry, CNumerics *numerics, CConfig *c } -void CFEASolver::BC_Clamped_Post(CGeometry *geometry, CNumerics *numerics, CConfig *config, unsigned short val_marker) { +void CFEASolver::BC_Clamped_Post(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker) { bool dynamic = config->GetTime_Domain(); - su2double zeros[3] = {0.0, 0.0, 0.0}; + su2double zeros[MAXNVAR] = {0.0}; for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { @@ -1791,7 +1798,70 @@ void CFEASolver::BC_Clamped_Post(CGeometry *geometry, CNumerics *numerics, CConf } -void CFEASolver::BC_DispDir(CGeometry *geometry, CNumerics *numerics, CConfig *config, unsigned short val_marker) { +void CFEASolver::BC_Sym_Plane(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker) { + + if (geometry->GetnElem_Bound(val_marker) == 0) return; + + const bool dynamic = config->GetTime_Domain(); + + /*--- Determine axis of symmetry based on the normal of the first element in the marker. ---*/ + + const su2double* nodeCoord[MAXNNODE_2D] = {nullptr}; + + const bool quad = (geometry->bound[val_marker][0]->GetVTK_Type() == QUADRILATERAL); + const unsigned short nNodes = quad? 4 : nDim; + + for (auto iNode = 0u; iNode < nNodes; iNode++) { + auto iPoint = geometry->bound[val_marker][0]->GetNode(iNode); + nodeCoord[iNode] = geometry->node[iPoint]->GetCoord(); + } + + su2double normal[MAXNDIM] = {0.0}; + + switch (nNodes) { + case 2: LineNormal(nodeCoord, normal); break; + case 3: TriangleNormal(nodeCoord, normal); break; + case 4: QuadrilateralNormal(nodeCoord, normal); break; + } + + auto axis = 0u; + for (auto iDim = 1u; iDim < MAXNDIM; ++iDim) + axis = (fabs(normal[iDim]) > fabs(normal[axis]))? iDim : axis; + + if (fabs(normal[axis]) < 0.99*Norm(int(MAXNDIM),normal)) { + SU2_MPI::Error("The structural solver only supports axis-aligned symmetry planes.",CURRENT_FUNCTION); + } + + /*--- Impose zero displacement perpendicular to the symmetry plane. ---*/ + + for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + /*--- Get node index ---*/ + const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode(); + + /*--- Set and enforce solution at current and previous time-step ---*/ + nodes->SetSolution(iPoint, axis, 0.0); + + if (dynamic) { + nodes->SetSolution_Vel(iPoint, axis, 0.0); + nodes->SetSolution_Accel(iPoint, axis, 0.0); + nodes->Set_Solution_time_n(iPoint, axis, 0.0); + nodes->SetSolution_Vel_time_n(iPoint, axis, 0.0); + nodes->SetSolution_Accel_time_n(iPoint, axis, 0.0); + } + + /*--- Set and enforce 0 solution for mesh deformation ---*/ + nodes->SetBound_Disp(iPoint, axis, 0.0); + + LinSysSol(iPoint, axis) = 0.0; + LinSysReact(iPoint, axis) = 0.0; + Jacobian.EnforceSolutionAtDOF(iPoint, axis, su2double(0.0), LinSysRes); + + } + +} + +void CFEASolver::BC_DispDir(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker) { unsigned short iDim; @@ -1799,13 +1869,7 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, CNumerics *numerics, CConfig *c su2double DispDirVal = config->GetDisp_Dir_Value(TagBound); su2double DispDirMult = config->GetDisp_Dir_Multiplier(TagBound); const su2double *DispDirLocal = config->GetDisp_Dir(TagBound); - - su2double DispDir[3] = {0.0}; - - su2double DispDirMod = 0.0; - for (iDim = 0; iDim < nDim; iDim++) - DispDirMod += DispDirLocal[iDim]*DispDirLocal[iDim]; - DispDirMod = sqrt(DispDirMod); + su2double DispDirMod = Norm(nDim, DispDirLocal); su2double CurrentTime = config->GetCurrent_DynTime(); su2double RampTime = config->GetRamp_Time(); @@ -1813,6 +1877,7 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, CNumerics *numerics, CConfig *c su2double TotalDisp = ModAmpl * DispDirVal * DispDirMult / DispDirMod; + su2double DispDir[MAXNVAR] = {0.0}; for (iDim = 0; iDim < nDim; iDim++) DispDir[iDim] = TotalDisp * DispDirLocal[iDim]; @@ -1835,81 +1900,118 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, CNumerics *numerics, CConfig *c void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, CNumerics **numerics, unsigned short iMesh) { - unsigned short iVar; - unsigned long iPoint, total_index; + const bool nonlinear_analysis = (config->GetGeometricConditions() == LARGE_DEFORMATIONS); - bool nonlinear_analysis = (config->GetGeometricConditions() == LARGE_DEFORMATIONS); - bool disc_adj_fem = (config->GetKind_Solver() == DISC_ADJ_FEM); + /*--- Compute stresses for monitoring and output. ---*/ - if (disc_adj_fem && nonlinear_analysis) { + Compute_NodalStress(geometry, numerics, config); - /*--- For nonlinear discrete adjoint, we have 3 convergence criteria ---*/ + /*--- Compute the objective function to be able to monitor it. ---*/ - /*--- UTOL = norm(Delta_U(k)): ABSOLUTE, norm of the incremental displacements ------------*/ - /*--- RTOL = norm(Residual(k): ABSOLUTE, norm of the residual (T-F) -----------------------*/ - /*--- ETOL = Delta_U(k) * Residual(k): ABSOLUTE, norm of the product disp * res -----------*/ + const auto kindObjFunc = config->GetKind_ObjFunc(); - Conv_Check[0] = LinSysSol.norm(); // Norm of the delta-solution vector - Conv_Check[1] = LinSysRes.norm(); // Norm of the residual - Conv_Check[2] = LinSysSol.dot(LinSysRes); // Position for the energy tolerance + if (((kindObjFunc == REFERENCE_GEOMETRY) || (kindObjFunc == REFERENCE_NODE)) && + ((config->GetDV_FEA() == YOUNG_MODULUS) || (config->GetDV_FEA() == DENSITY_VAL))) { + Stiffness_Penalty(geometry, solver_container, numerics, config); } - if (!nonlinear_analysis) { + switch (kindObjFunc) { + case REFERENCE_GEOMETRY: Compute_OFRefGeom(geometry, config); break; + case REFERENCE_NODE: Compute_OFRefNode(geometry, config); break; + case VOLUME_FRACTION: Compute_OFVolFrac(geometry, config); break; + case TOPOL_DISCRETENESS: Compute_OFVolFrac(geometry, config); break; + case TOPOL_COMPLIANCE: Compute_OFCompliance(geometry, config); break; + } + + if (nonlinear_analysis) { + + /*--- For nonlinear analysis we have 3 convergence criteria: ---*/ + /*--- UTOL = norm(Delta_U(k)): ABSOLUTE, norm of the incremental displacements ---*/ + /*--- RTOL = norm(Residual(k): ABSOLUTE, norm of the residual (T-F) ---*/ + /*--- ETOL = Delta_U(k) * Residual(k): ABSOLUTE, energy norm ---*/ + + SU2_OMP_PARALLEL + { + su2double utol = LinSysSol.norm(); + su2double rtol = LinSysRes.norm(); + su2double etol = LinSysSol.dot(LinSysRes); + + SU2_OMP_MASTER + { + Conv_Check[0] = utol; + Conv_Check[1] = rtol; + Conv_Check[2] = etol; + } + } // end parallel + } + else { /*--- If the problem is linear, the only check we do is the RMS of the residuals. ---*/ /*--- Compute the residual Ax-f ---*/ +#ifndef CODI_FORWARD_TYPE + CSysVector LinSysAux(nPoint, nPointDomain, nVar, nullptr); +#else + CSysVector LinSysAux(nPoint, nPointDomain, nVar, nullptr); +#endif + SU2_OMP_PARALLEL { #ifndef CODI_REVERSE_TYPE - Jacobian.ComputeResidual(LinSysSol, LinSysRes, LinSysAux); + Jacobian.ComputeResidual(LinSysSol, LinSysRes, LinSysAux); #else - /*--- We need temporaries to interface with the matrix ---*/ - { - CSysVector sol, res; - sol.PassiveCopy(LinSysSol); - res.PassiveCopy(LinSysRes); - CSysVector aux(res); - Jacobian.ComputeResidual(sol, res, aux); - LinSysAux.PassiveCopy(aux); - } + /*--- We need temporaries to interface with the passive matrix. ---*/ + CSysVector sol, res; + sol.PassiveCopy(LinSysSol); + res.PassiveCopy(LinSysRes); + Jacobian.ComputeResidual(sol, res, LinSysAux); #endif - } // end SU2_OMP_PARALLEL /*--- Set maximum residual to zero. ---*/ - for (iVar = 0; iVar < nVar; iVar++) { + SU2_OMP_MASTER + for (auto iVar = 0ul; iVar < nVar; iVar++) { SetRes_RMS(iVar, 0.0); SetRes_Max(iVar, 0.0, 0); } /*--- Compute the residual. ---*/ - for (iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { - for (iVar = 0; iVar < nVar; iVar++) { - total_index = iPoint*nVar+iVar; - AddRes_RMS(iVar, LinSysAux[total_index]*LinSysAux[total_index]); - AddRes_Max(iVar, fabs(LinSysAux[total_index]), - geometry->node[iPoint]->GetGlobalIndex(), - geometry->node[iPoint]->GetCoord()); + su2double resMax[MAXNVAR] = {0.0}, resRMS[MAXNVAR] = {0.0}; + const su2double* coordMax[MAXNVAR] = {nullptr}; + unsigned long idxMax[MAXNVAR] = {0}; + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { + for (auto iVar = 0ul; iVar < nVar; iVar++) { + su2double Res = fabs(LinSysAux(iPoint, iVar)); + resRMS[iVar] += Res*Res; + if (Res > resMax[iVar]) { + resMax[iVar] = Res; + idxMax[iVar] = iPoint; + coordMax[iVar] = geometry->node[iPoint]->GetCoord(); + } } } + SU2_OMP_CRITICAL + for (auto iVar = 0ul; iVar < nVar; iVar++) { + AddRes_RMS(iVar, resRMS[iVar]); + AddRes_Max(iVar, resMax[iVar], geometry->node[idxMax[iVar]]->GetGlobalIndex(), coordMax[iVar]); + } + SU2_OMP_BARRIER /*--- Compute the root mean square residual. ---*/ - + SU2_OMP_MASTER SetResidual_RMS(geometry, config); - } - - /*--- MPI solution (common to every mode). ---*/ + } // end SU2_OMP_PARALLEL - InitiateComms(geometry, config, SOLUTION_FEA); - CompleteComms(geometry, config, SOLUTION_FEA); + } } -void CFEASolver::BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, CConfig *config, unsigned short val_marker) { +void CFEASolver::BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker) { /*--- Determine whether the load conditions are applied in the reference or in the current configuration. ---*/ @@ -1932,13 +2034,9 @@ void CFEASolver::BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, CConfi for (unsigned long iElem = 0; iElem < geometry->GetnElem_Bound(val_marker); iElem++) { unsigned short iNode, iDim; - unsigned long indexNode[4] = {0,0,0,0}; - - su2double nodeCoord_ref[4][3], nodeCoord_curr[4][3]; - - for (iNode = 0; iNode < 4; iNode++) - for (iDim = 0; iDim < 3; iDim++) - nodeCoord_ref[iNode][iDim] = nodeCoord_curr[iNode][iDim] = 0.0; + unsigned long indexNode[MAXNNODE_2D] = {0}; + su2double nodeCoord_ref[MAXNNODE_2D][MAXNDIM] = {{0.0}}; + su2double nodeCoord_curr[MAXNNODE_2D][MAXNDIM] = {{0.0}}; /*--- Identify the kind of boundary element. ---*/ @@ -1964,8 +2062,8 @@ void CFEASolver::BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, CConfi /*--- Compute area vectors in reference and current configurations. ---*/ - su2double normal_ref[3] = {0.0, 0.0, 0.0}; - su2double normal_curr[3] = {0.0, 0.0, 0.0}; + su2double normal_ref[MAXNDIM] = {0.0}; + su2double normal_curr[MAXNDIM] = {0.0}; switch (nNodes) { case 2: LineNormal(nodeCoord_ref, normal_ref); break; @@ -2007,7 +2105,7 @@ void CFEASolver::BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, CConfi } -void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, CConfig *config, unsigned short val_marker) { +void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker) { auto TagBound = config->GetMarker_All_TagBound(val_marker); su2double LoadDirVal = config->GetLoad_Dir_Value(TagBound); @@ -2015,23 +2113,21 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, CConfig * const su2double* Load_Dir_Local = config->GetLoad_Dir(TagBound); /*--- Compute the norm of the vector that was passed in the config file. ---*/ - su2double Norm = pow(Load_Dir_Local[0],2) + pow(Load_Dir_Local[1],2); - if (nDim==3) Norm += pow(Load_Dir_Local[2],2); - Norm = sqrt(Norm); + su2double LoadNorm = Norm(nDim, Load_Dir_Local); su2double CurrentTime=config->GetCurrent_DynTime(); su2double Ramp_Time = config->GetRamp_Time(); su2double ModAmpl = Compute_LoadCoefficient(CurrentTime, Ramp_Time, config); - const su2double TotalLoad = ModAmpl * LoadDirVal * LoadDirMult / Norm; + const su2double TotalLoad = ModAmpl * LoadDirVal * LoadDirMult / LoadNorm; for (unsigned long iElem = 0; iElem < geometry->GetnElem_Bound(val_marker); iElem++) { unsigned short iNode, iDim; - unsigned long indexNode[4] = {0,0,0,0}; + unsigned long indexNode[MAXNNODE_2D] = {0}; - const su2double* nodeCoord[4] = {nullptr, nullptr, nullptr, nullptr}; + const su2double* nodeCoord[MAXNNODE_2D] = {nullptr}; /*--- Identify the kind of boundary element. ---*/ @@ -2047,7 +2143,7 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, CConfig * /*--- Compute area of the boundary element. ---*/ - su2double normal[3] = {0.0, 0.0, 0.0}; + su2double normal[MAXNDIM] = {0.0}; switch (nNodes) { case 2: LineNormal(nodeCoord, normal); break; @@ -2055,7 +2151,7 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, CConfig * case 4: QuadrilateralNormal(nodeCoord, normal); break; } - su2double area = sqrt(pow(normal[0],2) + pow(normal[1],2) + pow(normal[2],2)); + su2double area = Norm(int(MAXNDIM),normal); /*--- Compute load vector and update surface load for each node of the boundary element. ---*/ @@ -2070,16 +2166,16 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, CConfig * } -void CFEASolver::BC_Damper(CGeometry *geometry, CNumerics *numerics, CConfig *config, unsigned short val_marker) { +void CFEASolver::BC_Damper(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker) { const su2double dampConst = config->GetDamper_Constant(config->GetMarker_All_TagBound(val_marker)); for (auto iElem = 0ul; iElem < geometry->GetnElem_Bound(val_marker); iElem++) { unsigned short iNode, iDim; - unsigned long indexNode[4] = {0}; + unsigned long indexNode[MAXNNODE_2D] = {0}; - su2double nodeCoord[4][3] = {{0.0}}; + su2double nodeCoord[MAXNNODE_2D][MAXNDIM] = {{0.0}}; bool quad = (geometry->bound[val_marker][iElem]->GetVTK_Type() == QUADRILATERAL); unsigned short nNodes = quad? 4 : nDim; @@ -2097,7 +2193,7 @@ void CFEASolver::BC_Damper(CGeometry *geometry, CNumerics *numerics, CConfig *co /*--- Compute the area of the surface element. ---*/ - su2double normal[3] = {0.0}; + su2double normal[MAXNDIM] = {0.0}; switch (nNodes) { case 2: LineNormal(nodeCoord, normal); break; @@ -2105,7 +2201,7 @@ void CFEASolver::BC_Damper(CGeometry *geometry, CNumerics *numerics, CConfig *co case 4: QuadrilateralNormal(nodeCoord, normal); break; } - su2double area = sqrt(pow(normal[0],2) + pow(normal[1],2) + pow(normal[2],2)); + su2double area = Norm(int(MAXNDIM),normal); /*--- Compute damping forces. ---*/ @@ -2127,7 +2223,7 @@ void CFEASolver::BC_Damper(CGeometry *geometry, CNumerics *numerics, CConfig *co } -void CFEASolver::BC_Deforming(CGeometry *geometry, CNumerics *numerics, CConfig *config, unsigned short val_marker){ +void CFEASolver::BC_Deforming(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker){ for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { @@ -2135,7 +2231,7 @@ void CFEASolver::BC_Deforming(CGeometry *geometry, CNumerics *numerics, CConfig auto iNode = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Retrieve the boundary displacement ---*/ - su2double Disp[3] = {0.0, 0.0, 0.0}; + su2double Disp[MAXNVAR] = {0.0}; for (unsigned short iDim = 0; iDim < nDim; iDim++) Disp[iDim] = nodes->GetBound_Disp(iNode,iDim); @@ -2170,8 +2266,8 @@ void CFEASolver::Integrate_FSI_Loads(CGeometry *geometry, const CConfig *config) for (auto iElem = 0u; iElem < geometry->GetnElem_Bound(iMarker); ++iElem) { /*--- Define the boundary element. ---*/ - unsigned long nodeList[4] = {0}; - su2double coords[4][3] = {{0.0}}; + unsigned long nodeList[MAXNNODE_2D] = {0}; + su2double coords[MAXNNODE_2D][MAXNDIM] = {{0.0}}; bool quad = geometry->bound[iMarker][iElem]->GetVTK_Type() == QUADRILATERAL; auto nNode = quad? 4u : nDim; @@ -2183,7 +2279,7 @@ void CFEASolver::Integrate_FSI_Loads(CGeometry *geometry, const CConfig *config) } /*--- Compute the area contribution to each node. ---*/ - su2double normal[3] = {0.0}; + su2double normal[MAXNDIM] = {0.0}; switch (nNode) { case 2: LineNormal(coords, normal); break; @@ -2191,7 +2287,7 @@ void CFEASolver::Integrate_FSI_Loads(CGeometry *geometry, const CConfig *config) case 4: QuadrilateralNormal(coords, normal); break; } - su2double area = sqrt(pow(normal[0],2) + pow(normal[1],2) + pow(normal[2],2)) / nNode; + su2double area = Norm(int(MAXNDIM),normal) / nNode; /*--- Update area of nodes. ---*/ for (auto iNode = 0u; iNode < nNode; ++iNode) { @@ -2209,7 +2305,7 @@ void CFEASolver::Integrate_FSI_Loads(CGeometry *geometry, const CConfig *config) for (auto it = vertexArea.begin(); it != vertexArea.end(); ++it) { auto iPoint = it->first; su2double area = it->second; - su2double force[3] = {0.0}; + su2double force[MAXNDIM] = {0.0}; for (auto iDim = 0u; iDim < nDim; ++iDim) force[iDim] = nodes->Get_FlowTraction(iPoint,iDim)*area; nodes->Set_FlowTraction(iPoint, force); @@ -2217,7 +2313,7 @@ void CFEASolver::Integrate_FSI_Loads(CGeometry *geometry, const CConfig *config) } -su2double CFEASolver::Compute_LoadCoefficient(su2double CurrentTime, su2double RampTime, CConfig *config){ +su2double CFEASolver::Compute_LoadCoefficient(su2double CurrentTime, su2double RampTime, const CConfig *config){ su2double LoadCoeff = 1.0; @@ -2290,7 +2386,7 @@ su2double CFEASolver::Compute_LoadCoefficient(su2double CurrentTime, su2double R } -void CFEASolver::ImplicitNewmark_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { +void CFEASolver::ImplicitNewmark_Iteration(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { const bool first_iter = (config->GetInnerIter() == 0); const bool dynamic = (config->GetTime_Domain()); @@ -2350,7 +2446,7 @@ void CFEASolver::ImplicitNewmark_Iteration(CGeometry *geometry, CSolver **solver * */ if ((nonlinear_analysis && (newton_raphson || first_iter)) || linear_analysis) { - Jacobian.MatrixMatrixAddition(a_dt[0], MassMatrix); + Jacobian.MatrixMatrixAddition(SU2_TYPE::GetValue(a_dt[0]), MassMatrix); } /*--- Loop over all points, and set aux vector TimeRes_Aux = a0*U+a2*U'+a3*U'' ---*/ @@ -2366,8 +2462,7 @@ void CFEASolver::ImplicitNewmark_Iteration(CGeometry *geometry, CSolver **solver } /*--- Add M*TimeRes_Aux to the residual. ---*/ - - MassMatrix.MatrixVectorProduct(TimeRes_Aux, TimeRes, geometry, config); + Compute_MassRes(geometry, numerics, config); LinSysRes += TimeRes; } @@ -2375,7 +2470,7 @@ void CFEASolver::ImplicitNewmark_Iteration(CGeometry *geometry, CSolver **solver } -void CFEASolver::ImplicitNewmark_Update(CGeometry *geometry, CSolver **solver_container, CConfig *config) { +void CFEASolver::ImplicitNewmark_Update(CGeometry *geometry, CConfig *config) { const bool dynamic = (config->GetTime_Domain()); @@ -2429,7 +2524,7 @@ void CFEASolver::ImplicitNewmark_Update(CGeometry *geometry, CSolver **solver_co } -void CFEASolver::ImplicitNewmark_Relaxation(CGeometry *geometry, CSolver **solver_container, CConfig *config) { +void CFEASolver::ImplicitNewmark_Relaxation(CGeometry *geometry, CConfig *config) { const bool dynamic = (config->GetTime_Domain()); @@ -2493,7 +2588,7 @@ void CFEASolver::ImplicitNewmark_Relaxation(CGeometry *geometry, CSolver **solve } -void CFEASolver::GeneralizedAlpha_Iteration(CGeometry *geometry, CSolver **solver_container, CConfig *config) { +void CFEASolver::GeneralizedAlpha_Iteration(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { const bool first_iter = (config->GetInnerIter() == 0); const bool dynamic = (config->GetTime_Domain()); @@ -2549,7 +2644,7 @@ void CFEASolver::GeneralizedAlpha_Iteration(CGeometry *geometry, CSolver **solve /*--- See notes on logic in ImplicitNewmark_Iteration(). ---*/ if ((nonlinear_analysis && (newton_raphson || first_iter)) || linear_analysis) { - Jacobian.MatrixMatrixAddition(a_dt[0], MassMatrix); + Jacobian.MatrixMatrixAddition(SU2_TYPE::GetValue(a_dt[0]), MassMatrix); } /*--- Loop over all points, and set aux vector TimeRes_Aux = a0*U+a2*U'+a3*U'' ---*/ @@ -2565,8 +2660,7 @@ void CFEASolver::GeneralizedAlpha_Iteration(CGeometry *geometry, CSolver **solve } /*--- Add M*TimeRes_Aux to the residual. ---*/ - - MassMatrix.MatrixVectorProduct(TimeRes_Aux, TimeRes, geometry, config); + Compute_MassRes(geometry, numerics, config); LinSysRes += TimeRes; SU2_OMP_BARRIER @@ -2602,7 +2696,7 @@ void CFEASolver::GeneralizedAlpha_Iteration(CGeometry *geometry, CSolver **solve } -void CFEASolver::GeneralizedAlpha_UpdateDisp(CGeometry *geometry, CSolver **solver_container, CConfig *config) { +void CFEASolver::GeneralizedAlpha_UpdateDisp(CGeometry *geometry, CConfig *config) { /*--- Update displacement components of the solution. ---*/ @@ -2613,12 +2707,12 @@ void CFEASolver::GeneralizedAlpha_UpdateDisp(CGeometry *geometry, CSolver **solv /*--- Perform the MPI communication of the solution, displacements only. ---*/ - InitiateComms(geometry, config, SOLUTION_DISPONLY); - CompleteComms(geometry, config, SOLUTION_DISPONLY); + InitiateComms(geometry, config, SOLUTION); + CompleteComms(geometry, config, SOLUTION); } -void CFEASolver::GeneralizedAlpha_UpdateSolution(CGeometry *geometry, CSolver **solver_container, CConfig *config) { +void CFEASolver::GeneralizedAlpha_UpdateSolution(CGeometry *geometry, CConfig *config) { const su2double alpha_f = config->Get_Int_Coeffs(2); const su2double alpha_m = config->Get_Int_Coeffs(3); @@ -2677,7 +2771,7 @@ void CFEASolver::GeneralizedAlpha_UpdateSolution(CGeometry *geometry, CSolver ** } -void CFEASolver::GeneralizedAlpha_UpdateLoads(CGeometry *geometry, CSolver **solver_container, CConfig *config) { +void CFEASolver::GeneralizedAlpha_UpdateLoads(CGeometry *geometry, const CConfig *config) { /*--- Set the load conditions of the time step n+1 as the load conditions for time step n ---*/ nodes->Set_SurfaceLoad_Res_n(); @@ -2716,14 +2810,10 @@ void CFEASolver::Solve_System(CGeometry *geometry, CConfig *config) { } -void CFEASolver::PredictStruct_Displacement(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) { - - const unsigned short predOrder = fea_config->GetPredictorOrder(); - const su2double Delta_t = fea_config->GetDelta_DynTime(); +void CFEASolver::PredictStruct_Displacement(CGeometry *geometry, CConfig *config) { - auto fea_nodes = fea_solution[MESH_0][FEA_SOL]->GetNodes(); + const unsigned short predOrder = config->GetPredictorOrder(); + const su2double Delta_t = config->GetDelta_DynTime(); if(predOrder > 2 && rank == MASTER_NODE) cout << "Higher order predictor not implemented. Solving with order 0." << endl; @@ -2736,9 +2826,9 @@ void CFEASolver::PredictStruct_Displacement(CGeometry **fea_geometry, switch (predOrder) { case 1: { - const su2double* solDisp = fea_nodes->GetSolution(iPoint); - const su2double* solVel = fea_nodes->GetSolution_Vel(iPoint); - su2double* valPred = fea_nodes->GetSolution_Pred(iPoint); + const su2double* solDisp = nodes->GetSolution(iPoint); + const su2double* solVel = nodes->GetSolution_Vel(iPoint); + su2double* valPred = nodes->GetSolution_Pred(iPoint); for (iDim=0; iDim < nDim; iDim++) { valPred[iDim] = solDisp[iDim] + Delta_t*solVel[iDim]; @@ -2746,10 +2836,10 @@ void CFEASolver::PredictStruct_Displacement(CGeometry **fea_geometry, } break; case 2: { - const su2double* solDisp = fea_nodes->GetSolution(iPoint); - const su2double* solVel = fea_nodes->GetSolution_Vel(iPoint); - const su2double* solVel_tn = fea_nodes->GetSolution_Vel_time_n(iPoint); - su2double* valPred = fea_nodes->GetSolution_Pred(iPoint); + const su2double* solDisp = nodes->GetSolution(iPoint); + const su2double* solVel = nodes->GetSolution_Vel(iPoint); + const su2double* solVel_tn = nodes->GetSolution_Vel_time_n(iPoint); + su2double* valPred = nodes->GetSolution_Pred(iPoint); for (iDim=0; iDim < nDim; iDim++) { valPred[iDim] = solDisp[iDim] + 0.5*Delta_t*(3*solVel[iDim]-solVel_tn[iDim]); @@ -2757,16 +2847,18 @@ void CFEASolver::PredictStruct_Displacement(CGeometry **fea_geometry, } break; default: { - fea_nodes->SetSolution_Pred(iPoint); + nodes->SetSolution_Pred(iPoint); } break; } } + InitiateComms(geometry, config, SOLUTION_PRED); + CompleteComms(geometry, config, SOLUTION_PRED); + } -void CFEASolver::ComputeAitken_Coefficient(CGeometry **fea_geometry, CConfig *fea_config, - CSolver ***fea_solution, unsigned long iOuterIter) { +void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, CConfig *config, unsigned long iOuterIter) { unsigned long iPoint, iDim; su2double rbuf_numAitk = 0, sbuf_numAitk = 0; @@ -2776,11 +2868,11 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry **fea_geometry, CConfig *fe const su2double *dispCalc = nullptr; const su2double *dispPred_Old = nullptr; const su2double *dispCalc_Old = nullptr; - su2double deltaU[3] = {0.0, 0.0, 0.0}, deltaU_p1[3] = {0.0, 0.0, 0.0}; - su2double delta_deltaU[3] = {0.0, 0.0, 0.0}; + su2double deltaU[MAXNVAR] = {0.0}, deltaU_p1[MAXNVAR] = {0.0}; + su2double delta_deltaU[MAXNVAR] = {0.0}; su2double WAitkDyn_tn1, WAitkDyn_Max, WAitkDyn_Min, WAitkDyn; - unsigned short RelaxMethod_FSI = fea_config->GetRelaxation_Method_FSI(); + unsigned short RelaxMethod_FSI = config->GetRelaxation_Method_FSI(); /*--- Only when there is movement, and a dynamic coefficient is requested, it makes sense to compute the Aitken's coefficient ---*/ @@ -2791,7 +2883,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry **fea_geometry, CConfig *fe } else if (RelaxMethod_FSI == FIXED_PARAMETER) { - SetWAitken_Dyn(fea_config->GetAitkenStatRelax()); + SetWAitken_Dyn(config->GetAitkenStatRelax()); } else if (RelaxMethod_FSI == AITKEN_DYNAMIC) { @@ -2799,8 +2891,8 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry **fea_geometry, CConfig *fe if (iOuterIter == 0) { WAitkDyn_tn1 = GetWAitken_Dyn_tn1(); - WAitkDyn_Max = fea_config->GetAitkenDynMaxInit(); - WAitkDyn_Min = fea_config->GetAitkenDynMinInit(); + WAitkDyn_Max = config->GetAitkenDynMaxInit(); + WAitkDyn_Min = config->GetAitkenDynMinInit(); WAitkDyn = min(WAitkDyn_tn1, WAitkDyn_Max); WAitkDyn = max(WAitkDyn, WAitkDyn_Min); @@ -2812,10 +2904,10 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry **fea_geometry, CConfig *fe // To nPointDomain; we need to communicate the values for (iPoint = 0; iPoint < nPointDomain; iPoint++) { - dispPred = fea_solution[MESH_0][FEA_SOL]->GetNodes()->GetSolution_Pred(iPoint); - dispPred_Old = fea_solution[MESH_0][FEA_SOL]->GetNodes()->GetSolution_Pred_Old(iPoint); - dispCalc = fea_solution[MESH_0][FEA_SOL]->GetNodes()->GetSolution(iPoint); - dispCalc_Old = fea_solution[MESH_0][FEA_SOL]->GetNodes()->GetSolution_Old(iPoint); + dispPred = nodes->GetSolution_Pred(iPoint); + dispPred_Old = nodes->GetSolution_Pred_Old(iPoint); + dispCalc = nodes->GetSolution(iPoint); + dispCalc_Old = nodes->GetSolution_Old(iPoint); for (iDim = 0; iDim < nDim; iDim++) { @@ -2857,9 +2949,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry **fea_geometry, CConfig *fe } -void CFEASolver::SetAitken_Relaxation(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) { +void CFEASolver::SetAitken_Relaxation(CGeometry *geometry, CConfig *config) { const su2double WAitken = GetWAitken_Dyn(); @@ -2868,14 +2958,14 @@ void CFEASolver::SetAitken_Relaxation(CGeometry **fea_geometry, for (unsigned long iPoint=0; iPoint < nPointDomain; iPoint++) { /*--- Retrieve pointers to the predicted and calculated solutions ---*/ - su2double* dispPred = fea_solution[MESH_0][FEA_SOL]->GetNodes()->GetSolution_Pred(iPoint); - const su2double* dispCalc = fea_solution[MESH_0][FEA_SOL]->GetNodes()->GetSolution(iPoint); + su2double* dispPred = nodes->GetSolution_Pred(iPoint); + const su2double* dispCalc = nodes->GetSolution(iPoint); /*--- Set predicted solution as the old predicted solution ---*/ - fea_solution[MESH_0][FEA_SOL]->GetNodes()->SetSolution_Pred_Old(iPoint); + nodes->SetSolution_Pred_Old(iPoint); /*--- Set calculated solution as the old solution (needed for dynamic Aitken relaxation) ---*/ - fea_solution[MESH_0][FEA_SOL]->GetNodes()->SetSolution_Old(iPoint, dispCalc); + nodes->SetSolution_Old(iPoint, dispCalc); /*--- Apply the Aitken relaxation ---*/ for (unsigned short iDim=0; iDim < nDim; iDim++) { @@ -2883,28 +2973,12 @@ void CFEASolver::SetAitken_Relaxation(CGeometry **fea_geometry, } } -} - -void CFEASolver::Update_StructSolution(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) { - - SU2_OMP_PARALLEL_(for schedule(static,omp_chunk_size)) - for (unsigned long iPoint=0; iPoint < nPointDomain; iPoint++) { - - auto valSolutionPred = fea_solution[MESH_0][FEA_SOL]->GetNodes()->GetSolution_Pred(iPoint); - - fea_solution[MESH_0][FEA_SOL]->GetNodes()->SetSolution(iPoint, valSolutionPred); - } - - /*--- Perform the MPI communication of the solution, displacements only ---*/ - - InitiateComms(fea_geometry[MESH_0], fea_config, SOLUTION_DISPONLY); - CompleteComms(fea_geometry[MESH_0], fea_config, SOLUTION_DISPONLY); + InitiateComms(geometry, config, SOLUTION_PRED_OLD); + CompleteComms(geometry, config, SOLUTION_PRED_OLD); } -void CFEASolver::OutputForwardModeGradient(CConfig *config, bool newFile, +void CFEASolver::OutputForwardModeGradient(const CConfig *config, bool newFile, su2double fun, su2double fun_avg, su2double der, su2double der_avg) const { if (rank != MASTER_NODE) return; @@ -2974,10 +3048,9 @@ void CFEASolver::OutputForwardModeGradient(CConfig *config, bool newFile, } -void CFEASolver::Compute_OFRefGeom(CGeometry *geometry, CSolver **solver_container, CConfig *config){ +void CFEASolver::Compute_OFRefGeom(CGeometry *geometry, const CConfig *config){ bool fsi = config->GetFSI_Simulation(); - bool dynamic = config->GetTime_Domain(); unsigned long TimeIter = config->GetTimeIter(); su2double objective_function = 0.0; @@ -3009,17 +3082,18 @@ void CFEASolver::Compute_OFRefGeom(CGeometry *geometry, CSolver **solver_contain Total_OFRefGeom += PenaltyValue; Global_OFRefGeom += Total_OFRefGeom; - su2double objective_function_averaged = Global_OFRefGeom / (TimeIter + 1.0 + EPS); - /// TODO: Temporary output files for the objective function. Will be integrated in the output once it is refactored. + /*--- To be accessible from the output. ---*/ + Total_OFCombo = Total_OFRefGeom; - if (rank != MASTER_NODE) return; + /// TODO: Temporary output files for the direct mode. - if (config->GetDirectDiff() != NO_DERIVATIVE) { + if ((rank == MASTER_NODE) && (config->GetDirectDiff() != NO_DERIVATIVE)) { /*--- Forward mode AD results. ---*/ su2double local_forward_gradient = SU2_TYPE::GetDerivative(Total_OFRefGeom); + su2double objective_function_averaged = Global_OFRefGeom / (TimeIter + 1.0 + EPS); if (fsi) Total_ForwardGradient = local_forward_gradient; else Total_ForwardGradient += local_forward_gradient; @@ -3029,32 +3103,15 @@ void CFEASolver::Compute_OFRefGeom(CGeometry *geometry, CSolver **solver_contain OutputForwardModeGradient(config, false, Total_OFRefGeom, objective_function_averaged, local_forward_gradient, averaged_gradient); } - else { - cout << "Objective function: " << Total_OFRefGeom << "." << endl; - ofstream myfile_res; - myfile_res.open ("of_refgeom.dat"); - myfile_res.precision(15); - if (dynamic) myfile_res << scientific << objective_function_averaged << endl; - else myfile_res << scientific << Total_OFRefGeom << endl; - myfile_res.close(); - if (fsi) { - ofstream myfile_his; - myfile_his.open ("history_refgeom.dat",ios::app); - myfile_his.precision(15); - myfile_his << scientific << Total_OFRefGeom << endl; - myfile_his.close(); - } - } } -void CFEASolver::Compute_OFRefNode(CGeometry *geometry, CSolver **solver_container, CConfig *config){ +void CFEASolver::Compute_OFRefNode(CGeometry *geometry, const CConfig *config){ bool fsi = config->GetFSI_Simulation(); - bool dynamic = config->GetTime_Domain(); unsigned long TimeIter = config->GetTimeIter(); - su2double dist[3] = {0.0, 0.0, 0.0}, dist_reduce[3]; + su2double dist[MAXNVAR] = {0.0}, dist_reduce[MAXNVAR]; /*--- Convert global point index to local. ---*/ long iPoint = geometry->GetGlobal_to_Local_Point(config->GetRefNode_ID()); @@ -3066,26 +3123,23 @@ void CFEASolver::Compute_OFRefNode(CGeometry *geometry, CSolver **solver_contain } } - SU2_MPI::Allreduce(dist, dist_reduce, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + SU2_MPI::Allreduce(dist, dist_reduce, MAXNVAR, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - Total_OFRefNode = 0.0; - for (unsigned short iVar = 0; iVar < 3; ++iVar) - Total_OFRefNode += pow(dist_reduce[iVar],2); - - Total_OFRefNode = config->GetRefNode_Penalty()*sqrt(Total_OFRefNode) + PenaltyValue; + Total_OFRefNode = config->GetRefNode_Penalty() * Norm(int(MAXNVAR),dist_reduce) + PenaltyValue; Global_OFRefNode += Total_OFRefNode; - su2double objective_function_averaged = Global_OFRefNode / (TimeIter + 1.0 + EPS); - /// TODO: Temporary output files for the objective function. Will be integrated in the output once it is refactored. + /*--- To be accessible from the output. ---*/ + Total_OFCombo = Total_OFRefNode; - if (rank != MASTER_NODE) return; + /// TODO: Temporary output files for the direct mode. - if (config->GetDirectDiff() != NO_DERIVATIVE) { + if ((rank == MASTER_NODE) && (config->GetDirectDiff() != NO_DERIVATIVE)) { /*--- Forward mode AD results. ---*/ su2double local_forward_gradient = SU2_TYPE::GetDerivative(Total_OFRefNode); + su2double objective_function_averaged = Global_OFRefNode / (TimeIter + 1.0 + EPS); if (fsi) Total_ForwardGradient = local_forward_gradient; else Total_ForwardGradient += local_forward_gradient; @@ -3095,30 +3149,10 @@ void CFEASolver::Compute_OFRefNode(CGeometry *geometry, CSolver **solver_contain OutputForwardModeGradient(config, false, Total_OFRefNode, objective_function_averaged, local_forward_gradient, averaged_gradient); } - else { - cout << "Objective function: " << Total_OFRefNode << "." << endl; - ofstream myfile_res; - myfile_res.open ("of_refnode.dat"); - myfile_res.precision(15); - if (dynamic) myfile_res << scientific << objective_function_averaged << endl; - else myfile_res << scientific << Total_OFRefNode << endl; - myfile_res.close(); - - ofstream myfile_his; - myfile_his.open ("history_refnode.dat",ios::app); - myfile_his.precision(15); - myfile_his << TimeIter << "\t"; - myfile_his << scientific << Total_OFRefNode << "\t"; - myfile_his << scientific << objective_function_averaged << "\t"; - myfile_his << scientific << dist_reduce[0] << "\t"; - myfile_his << scientific << dist_reduce[1] << "\t"; - myfile_his << scientific << dist_reduce[2] << endl; - myfile_his.close(); - } } -void CFEASolver::Compute_OFVolFrac(CGeometry *geometry, CSolver **solver_container, CConfig *config) +void CFEASolver::Compute_OFVolFrac(CGeometry *geometry, const CConfig *config) { /*--- Perform a volume average of the physical density of the elements for topology optimization ---*/ @@ -3157,24 +3191,12 @@ void CFEASolver::Compute_OFVolFrac(CGeometry *geometry, CSolver **solver_contain else Total_OFVolFrac = integral/total_volume; - /// TODO: Temporary output files for the objective function. Will be integrated in the output once it is refactored. - - if (rank == MASTER_NODE) { - cout << "Objective function: " << Total_OFVolFrac << "." << endl; - - ofstream myfile_res; - if (config->GetKind_ObjFunc() == TOPOL_DISCRETENESS) - myfile_res.open ("of_topdisc.dat"); - else - myfile_res.open ("of_volfrac.dat"); + /*--- To be accessible from the output. ---*/ + Total_OFCombo = Total_OFVolFrac; - myfile_res.precision(15); - myfile_res << scientific << Total_OFVolFrac << endl; - myfile_res.close(); - } } -void CFEASolver::Compute_OFCompliance(CGeometry *geometry, CSolver **solver_container, CConfig *config) +void CFEASolver::Compute_OFCompliance(CGeometry *geometry, const CConfig *config) { /*--- Types of loads to consider ---*/ const bool body_forces = config->GetDeadLoad(); @@ -3224,21 +3246,19 @@ void CFEASolver::Compute_OFCompliance(CGeometry *geometry, CSolver **solver_cont SU2_MPI::Allreduce(&compliance, &Total_OFCompliance, 1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); - /// TODO: Temporary output files for the objective function. Will be integrated in the output once it is refactored. + /*--- To be accessible from the output. ---*/ + Total_OFCombo = Total_OFCompliance; - if (rank == MASTER_NODE) { - cout << "Objective function: " << Total_OFCompliance << "." << endl; - - ofstream file; - file.open("of_topcomp.dat"); - file.precision(15); - file << scientific << Total_OFCompliance << endl; - file.close(); - } } void CFEASolver::Stiffness_Penalty(CGeometry *geometry, CSolver **solver, CNumerics **numerics, CConfig *config){ + if (config->GetTotalDV_Penalty() == 0.0) { + /*--- No need to go into expensive computations. ---*/ + PenaltyValue = 0.0; + return; + } + su2double weightedValue = 0.0; su2double weightedValue_reduce = 0.0; su2double totalVolume = 0.0; @@ -3256,7 +3276,7 @@ void CFEASolver::Stiffness_Penalty(CGeometry *geometry, CSolver **solver, CNumer int EL_KIND; unsigned short iNode, nNodes, iDim; - unsigned long indexNode[MAXNNODE]; + unsigned long indexNode[MAXNNODE_3D]; GetElemKindAndNumNodes(geometry->elem[iElem]->GetVTK_Type(), EL_KIND, nNodes); @@ -3460,7 +3480,7 @@ void CFEASolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config) } -void CFEASolver::FilterElementDensities(CGeometry *geometry, CConfig *config) +void CFEASolver::FilterElementDensities(CGeometry *geometry, const CConfig *config) { /*--- Apply a filter to the design densities of the elements to generate the physical densities which are the ones used to penalize their stiffness. ---*/ diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 8fcd4392e01c..517f3dcf1e40 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -284,11 +284,9 @@ void CMeshSolver::SetMinMaxVolume(CGeometry *geometry, CConfig *config, bool upd void CMeshSolver::SetWallDistance(CGeometry *geometry, CConfig *config) { - unsigned long nVertex_SolidWall, ii, jj, iVertex, iPoint, pointID, iElem; - unsigned short iNodes, nNodes = 0; + unsigned long nVertex_SolidWall, ii, jj, iVertex, iPoint, pointID; unsigned short iMarker, iDim; su2double dist, MaxDistance_Local, MinDistance_Local; - su2double nodeDist, ElemDist; int rankID; /*--- Initialize min and max distance ---*/ @@ -345,8 +343,7 @@ void CMeshSolver::SetWallDistance(CGeometry *geometry, CConfig *config) { /*--- No solid wall boundary nodes in the entire mesh. Set the wall distance to zero for all nodes. ---*/ - for (iPoint=0; iPointGetnPoint(); ++iPoint) - geometry->node[iPoint]->SetWall_Distance(0.0); + geometry->SetWallDistance(0.0); } else { @@ -375,29 +372,34 @@ void CMeshSolver::SetWallDistance(CGeometry *geometry, CConfig *config) { } + SU2_OMP_PARALLEL + { /*--- Normalize distance from 0 to 1 ---*/ - for (iPoint=0; iPoint < nPoint; ++iPoint) { - nodeDist = nodes->GetWallDistance(iPoint)/MaxDistance; + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { + su2double nodeDist = nodes->GetWallDistance(iPoint)/MaxDistance; nodes->SetWallDistance(iPoint,nodeDist); } /*--- Compute the element distances ---*/ - for (iElem = 0; iElem < nElement; iElem++) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iElem = 0ul; iElem < nElement; iElem++) { int EL_KIND; + unsigned short nNodes = 0; GetElemKindAndNumNodes(geometry->elem[iElem]->GetVTK_Type(), EL_KIND, nNodes); /*--- Average the distance of the nodes in the element ---*/ - ElemDist = 0.0; - for (iNodes = 0; iNodes < nNodes; iNodes++){ - iPoint = geometry->elem[iElem]->GetNode(iNodes); + su2double ElemDist = 0.0; + for (auto iNode = 0u; iNode < nNodes; iNode++) { + auto iPoint = geometry->elem[iElem]->GetNode(iNode); ElemDist += nodes->GetWallDistance(iPoint); } ElemDist = ElemDist/su2double(nNodes); element[iElem].SetWallDistance(ElemDist); - + } } } @@ -461,29 +463,31 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig if (multizone) nodes->Set_BGSSolution_k(); - /*--- Compute the stiffness matrix. ---*/ - Compute_StiffMatrix(geometry[MESH_0], numerics, config); - - /*--- Clean residual, we do not want an incremental solution. ---*/ - SU2_OMP_PARALLEL - { - LinSysRes.SetValZero(); - } - - /*--- LinSysSol contains the non-transformed displacements in the periodic halo cells. - Hence we still need a communication of the transformed coordinates, otherwise periodicity - is not maintained. ---*/ + /*--- Capture a few MPI dependencies for AD. ---*/ geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, COORDINATES); geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, COORDINATES); - /*--- In the same way, communicate the displacements in the solver to make sure the halo - nodes receive the correct value of the displacement. ---*/ InitiateComms(geometry[MESH_0], config, SOLUTION); CompleteComms(geometry[MESH_0], config, SOLUTION); InitiateComms(geometry[MESH_0], config, MESH_DISPLACEMENTS); CompleteComms(geometry[MESH_0], config, MESH_DISPLACEMENTS); + /*--- Compute the stiffness matrix, no point recording because we clear the residual. ---*/ + + const bool ActiveTape = AD::TapeActive(); + AD::StopRecording(); + + Compute_StiffMatrix(geometry[MESH_0], numerics, config); + + if (ActiveTape) AD::StartRecording(); + + /*--- Clear residual (loses AD info), we do not want an incremental solution. ---*/ + SU2_OMP_PARALLEL + { + LinSysRes.SetValZero(); + } + /*--- Impose boundary conditions (all of them are ESSENTIAL BC's - displacements). ---*/ SetBoundaryDisplacements(geometry[MESH_0], numerics[FEA_TERM], config); @@ -498,8 +502,10 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig UpdateDualGrid(geometry[MESH_0], config); /*--- Check for failed deformation (negative volumes). ---*/ - /*--- In order to do this, we recompute the minimum and maximum area/volume for the mesh using the current coordinates. ---*/ + /*--- This is not recorded as it does not influence the solution. ---*/ + AD::StopRecording(); SetMinMaxVolume(geometry[MESH_0], config, true); + if (ActiveTape) AD::StartRecording(); /*--- The Grid Velocity is only computed if the problem is time domain ---*/ if (time_domain) ComputeGridVelocity(geometry[MESH_0], config); @@ -515,11 +521,10 @@ void CMeshSolver::UpdateGridCoord(CGeometry *geometry, CConfig *config){ /*--- LinSysSol contains the absolute x, y, z displacements. ---*/ SU2_OMP_PARALLEL_(for schedule(static,omp_chunk_size)) - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){ + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++){ for (unsigned short iDim = 0; iDim < nDim; iDim++) { - auto total_index = iPoint*nDim + iDim; /*--- Retrieve the displacement from the solution of the linear system ---*/ - su2double val_disp = LinSysSol[total_index]; + su2double val_disp = LinSysSol(iPoint, iDim); /*--- Store the displacement of the mesh node ---*/ nodes->SetSolution(iPoint, iDim, val_disp); /*--- Compute the current coordinate as Mesh_Coord + Displacement ---*/ @@ -529,14 +534,10 @@ void CMeshSolver::UpdateGridCoord(CGeometry *geometry, CConfig *config){ } } - /*--- LinSysSol contains the non-transformed displacements in the periodic halo cells. - Hence we still need a communication of the transformed coordinates, otherwise periodicity - is not maintained. ---*/ + /*--- Communicate the updated displacements and mesh coordinates. ---*/ geometry->InitiateComms(geometry, config, COORDINATES); geometry->CompleteComms(geometry, config, COORDINATES); - /*--- In the same way, communicate the displacements in the solver to make sure the halo - nodes receive the correct value of the displacement. ---*/ InitiateComms(geometry, config, SOLUTION); CompleteComms(geometry, config, SOLUTION); @@ -668,6 +669,22 @@ void CMeshSolver::SetBoundaryDisplacements(CGeometry *geometry, CNumerics *numer } } + /*--- Clamp far away nodes according to deform limit. ---*/ + if ((config->GetDeform_Stiffness_Type() == SOLID_WALL_DISTANCE) && + (config->GetDeform_Limit() < MaxDistance)) { + + const su2double limit = config->GetDeform_Limit() / MaxDistance; + + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { + if (nodes->GetWallDistance(iPoint) <= limit) continue; + + su2double zeros[MAXNVAR] = {0.0}; + nodes->SetSolution(iPoint, zeros); + LinSysSol.SetBlock(iPoint, zeros); + Jacobian.EnforceSolutionAtNode(iPoint, zeros, LinSysRes); + } + } + } void CMeshSolver::SetDualTime_Mesh(void){ diff --git a/SU2_CFD/src/solvers/CRadP1Solver.cpp b/SU2_CFD/src/solvers/CRadP1Solver.cpp index acd84e2589ce..b2137ede4e0d 100644 --- a/SU2_CFD/src/solvers/CRadP1Solver.cpp +++ b/SU2_CFD/src/solvers/CRadP1Solver.cpp @@ -82,7 +82,6 @@ CRadP1Solver::CRadP1Solver(CGeometry* geometry, CConfig *config) : CRadSolver(ge LinSysSol.Initialize(nPoint, nPointDomain, nVar, 0.0); LinSysRes.Initialize(nPoint, nPointDomain, nVar, 0.0); - LinSysAux.Initialize(nPoint, nPointDomain, nVar, 0.0); /*--- Read farfield conditions from config ---*/ Temperature_Inf = config->GetTemperature_FreeStreamND(); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 92f9b81a4e60..f6c1ee273446 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1958,10 +1958,6 @@ void CSolver::InitiateComms(CGeometry *geometry, COUNT_PER_POINT = nVar*3; MPI_TYPE = COMM_TYPE_DOUBLE; break; - case SOLUTION_DISPONLY: - COUNT_PER_POINT = nVar; - MPI_TYPE = COMM_TYPE_DOUBLE; - break; case SOLUTION_PRED: COUNT_PER_POINT = nVar; MPI_TYPE = COMM_TYPE_DOUBLE; @@ -2103,10 +2099,6 @@ void CSolver::InitiateComms(CGeometry *geometry, bufDSend[buf_offset+nVar*2+iVar] = base_nodes->GetSolution_Accel_time_n(iPoint, iVar); } break; - case SOLUTION_DISPONLY: - for (iVar = 0; iVar < nVar; iVar++) - bufDSend[buf_offset+iVar] = base_nodes->GetSolution(iPoint, iVar); - break; case SOLUTION_PRED: for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution_Pred(iPoint, iVar); @@ -2271,10 +2263,6 @@ void CSolver::CompleteComms(CGeometry *geometry, base_nodes->SetSolution_Accel_time_n(iPoint, iVar, bufDRecv[buf_offset+nVar*2+iVar]); } break; - case SOLUTION_DISPONLY: - for (iVar = 0; iVar < nVar; iVar++) - base_nodes->SetSolution(iPoint, iVar, bufDRecv[buf_offset+iVar]); - break; case SOLUTION_PRED: for (iVar = 0; iVar < nVar; iVar++) base_nodes->SetSolution_Pred(iPoint, iVar, bufDRecv[buf_offset+iVar]); diff --git a/SU2_PY/topology_optimization.py b/SU2_PY/topology_optimization.py index 8131dd8b1980..8740fefa95b9 100755 --- a/SU2_PY/topology_optimization.py +++ b/SU2_PY/topology_optimization.py @@ -45,9 +45,9 @@ ####### SETUP ####### -obj_scale = 1/0.0125 # scale the objective so that it starts at 1-4 -con_scale = 1/0.5 # 1 over upper bound (e.g. max volume) -var_scale = 1.0 # variable scale +obj_scale = 1/1.25e-3 # scale the objective so that it starts at 1-4 +con_scale = 1/0.5 # 1 over upper bound (e.g. max volume) +var_scale = 1.0 # variable scale # maximum number of iterations maxJev_t = 1000 @@ -73,11 +73,10 @@ inputFile = "element_properties.dat" # names of the output files [objective value, objective gradient, constraint value, ...] -outputFiles = ["of_refnode.dat", "grad_ref_node.dat", - "of_volfrac.dat", "grad_vol_frac.dat"] +outputFiles = ["grad_compliance.dat", "grad_vol_frac.dat"] # settings for direct run and adjoint of the objective and constraint -fnames = ["settings.cfg", "settings_refnode.cfg", "settings_volfrac.cfg"] +fnames = ["settings.cfg", "settings_compliance.cfg", "settings_volfrac.cfg"] # use the DILATE, ERODE, (DILATE,ERODE), or (ERODE,DILATE) filters, the first value is # for gray initialization, then it is ramped until a solid-void topology is obtained @@ -89,10 +88,10 @@ class Driver: def __init__(self,commands,inputFile,configFiles,outputFiles): self._inputFile = inputFile - self._objValFile = outputFiles[0] - self._objDerFile = outputFiles[1] - self._conValFile = outputFiles[2] - self._conDerFile = outputFiles[3] + self._objValFile = "history.csv" + self._objDerFile = outputFiles[0] + self._conValFile = "history.csv" + self._conDerFile = outputFiles[1] self._objValCommand = commands[0]+configFiles[0]+" > objval.stdout" self._objDerCommand = commands[1]+configFiles[1]+" > objder.stdout" self._conDerCommand = commands[1]+configFiles[2]+" > conval.stdout" @@ -123,7 +122,7 @@ def obj_val(self,x): try: sp.call(self._objValCommand,shell=True) - fid = open(self._objValFile,"r"); val = float(fid.read()); fid.close() + fid = open(self._objValFile,"r"); val = float(fid.readlines()[1]); fid.close() # the return code of mpirun is useless, we test the value of the function self._assert_isfinite(val) except: @@ -170,7 +169,7 @@ def con_val(self,x): try: sp.call(self._conDerCommand,shell=True) - fid = open(self._conValFile,"r"); val = float(fid.read()); fid.close() + fid = open(self._conValFile,"r"); val = float(fid.readlines()[1]); fid.close() self._assert_isfinite(val) except: raise RuntimeError("Constraint function evaluation failed") diff --git a/TestCases/fea_fsi/MixElemsKnowles/config.cfg b/TestCases/fea_fsi/MixElemsKnowles/config.cfg index 90d8591abd91..fdb5428667b2 100644 --- a/TestCases/fea_fsi/MixElemsKnowles/config.cfg +++ b/TestCases/fea_fsi/MixElemsKnowles/config.cfg @@ -51,4 +51,4 @@ TABULAR_FORMAT= CSV RESTART_SOL= YES SOLUTION_FILENAME= solution_structure.dat -SCREEN_OUTPUT=(INNER_ITER, RMS_UTOL, RMS_RTOL, RMS_ETOL, VMS) +SCREEN_OUTPUT=(INNER_ITER, RMS_UTOL, RMS_RTOL, RMS_ETOL, COMBO, VMS) diff --git a/TestCases/fea_topology/quick_start/settings.cfg b/TestCases/fea_topology/quick_start/settings.cfg index 4c0c7c89b3b6..e1ce988d48d3 100644 --- a/TestCases/fea_topology/quick_start/settings.cfg +++ b/TestCases/fea_topology/quick_start/settings.cfg @@ -10,17 +10,13 @@ TOPOLOGY_OPTIMIZATION= YES FEA_FILENAME= element_properties.dat TOPOL_OPTIM_FILTER_RADIUS= 0.013 TOPOL_OPTIM_FILTER_KERNEL= ( DILATE,ERODE ) -TOPOL_OPTIM_KERNEL_PARAM= 200 +TOPOL_OPTIM_KERNEL_PARAM= 0.01 TOPOL_OPTIM_PROJECTION_TYPE= NO_PROJECTION TOPOL_OPTIM_PROJECTION_PARAM= 0.0 TOPOL_OPTIM_SIMP_EXPONENT= 3.0 TOPOL_OPTIM_OUTFILE= grad_ref_node.dat % -OBJECTIVE_FUNCTION= REFERENCE_NODE -REFERENCE_NODE= 890 -REFERENCE_NODE_DISPLACEMENT= (0.0, 0.0) -REFERENCE_NODE_PENALTY= 1.0 -DESIGN_VARIABLE_FEA= YOUNG_MODULUS +OBJECTIVE_FUNCTION= TOPOL_COMPLIANCE % % Numerics GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS @@ -47,13 +43,7 @@ LINEAR_SOLVER_ITER= 1000 MESH_FILENAME= mesh.su2 MESH_FORMAT= SU2 TABULAR_FORMAT= CSV -CONV_FILENAME= history -VOLUME_STRUCTURE_FILENAME= direct -VOLUME_ADJ_STRUCTURE_FILENAME= adjoint -SOLUTION_STRUCTURE_FILENAME= direct.dat -RESTART_STRUCTURE_FILENAME= direct.dat -WRT_SOL_FREQ= 9999 -WRT_CON_FREQ= 1 -WRT_CON_FREQ_DUALTIME= 1 -RESTART_SOL= NO - +SOLUTION_FILENAME= direct.dat +RESTART_FILENAME= direct.dat +OUTPUT_FILES= RESTART, PARAVIEW +HISTORY_OUTPUT= COMBO diff --git a/TestCases/fea_topology/quick_start/settings_refnode.cfg b/TestCases/fea_topology/quick_start/settings_compliance.cfg similarity index 66% rename from TestCases/fea_topology/quick_start/settings_refnode.cfg rename to TestCases/fea_topology/quick_start/settings_compliance.cfg index 90247d60e9e4..6c7df4e5f55b 100644 --- a/TestCases/fea_topology/quick_start/settings_refnode.cfg +++ b/TestCases/fea_topology/quick_start/settings_compliance.cfg @@ -10,17 +10,13 @@ TOPOLOGY_OPTIMIZATION= YES FEA_FILENAME= element_properties.dat TOPOL_OPTIM_FILTER_RADIUS= 0.013 TOPOL_OPTIM_FILTER_KERNEL= ( DILATE,ERODE ) -TOPOL_OPTIM_KERNEL_PARAM= 200 +TOPOL_OPTIM_KERNEL_PARAM= 0.01 TOPOL_OPTIM_PROJECTION_TYPE= NO_PROJECTION TOPOL_OPTIM_PROJECTION_PARAM= 0.0 TOPOL_OPTIM_SIMP_EXPONENT= 3.0 -TOPOL_OPTIM_OUTFILE= grad_ref_node.dat +TOPOL_OPTIM_OUTFILE= grad_compliance.dat % -OBJECTIVE_FUNCTION= REFERENCE_NODE -REFERENCE_NODE= 890 -REFERENCE_NODE_DISPLACEMENT= (0.0, 0.0) -REFERENCE_NODE_PENALTY= 1.0 -DESIGN_VARIABLE_FEA= YOUNG_MODULUS +OBJECTIVE_FUNCTION= TOPOL_COMPLIANCE % % Numerics GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS @@ -49,13 +45,8 @@ DISCADJ_LIN_PREC= ILU MESH_FILENAME= mesh.su2 MESH_FORMAT= SU2 TABULAR_FORMAT= CSV -CONV_FILENAME= history -VOLUME_STRUCTURE_FILENAME= direct -VOLUME_ADJ_STRUCTURE_FILENAME= adjoint -SOLUTION_STRUCTURE_FILENAME= direct.dat -RESTART_STRUCTURE_FILENAME= direct.dat -WRT_SOL_FREQ= 9999 -WRT_CON_FREQ= 1 -WRT_CON_FREQ_DUALTIME= 1 -RESTART_SOL= NO +SOLUTION_FILENAME= direct.dat +RESTART_FILENAME= direct.dat +OUTPUT_FILES= PARAVIEW +HISTORY_OUTPUT= COMBO diff --git a/TestCases/fea_topology/quick_start/settings_volfrac.cfg b/TestCases/fea_topology/quick_start/settings_volfrac.cfg index 1e5848fabe08..96dc122c3f66 100644 --- a/TestCases/fea_topology/quick_start/settings_volfrac.cfg +++ b/TestCases/fea_topology/quick_start/settings_volfrac.cfg @@ -10,14 +10,13 @@ TOPOLOGY_OPTIMIZATION= YES FEA_FILENAME= element_properties.dat TOPOL_OPTIM_FILTER_RADIUS= 0.013 TOPOL_OPTIM_FILTER_KERNEL= ( DILATE,ERODE ) -TOPOL_OPTIM_KERNEL_PARAM= 200 +TOPOL_OPTIM_KERNEL_PARAM= 0.01 TOPOL_OPTIM_PROJECTION_TYPE= NO_PROJECTION TOPOL_OPTIM_PROJECTION_PARAM= 0.0 TOPOL_OPTIM_SIMP_EXPONENT= 3.0 TOPOL_OPTIM_OUTFILE= grad_vol_frac.dat % OBJECTIVE_FUNCTION= VOLUME_FRACTION -DESIGN_VARIABLE_FEA= YOUNG_MODULUS % % Numerics GEOMETRIC_CONDITIONS= SMALL_DEFORMATIONS @@ -36,23 +35,18 @@ MARKER_LOAD= ( load, 1, 1, 0, -1, 0) % % Linear solver LINEAR_SOLVER= CONJUGATE_GRADIENT -LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_PREC= JACOBI LINEAR_SOLVER_ERROR= 1E-9 -LINEAR_SOLVER_ITER= 1000 +LINEAR_SOLVER_ITER= 1 DISCADJ_LIN_SOLVER= CONJUGATE_GRADIENT -DISCADJ_LIN_PREC= ILU +DISCADJ_LIN_PREC= JACOBI % % In/Out MESH_FILENAME= mesh.su2 MESH_FORMAT= SU2 TABULAR_FORMAT= CSV -CONV_FILENAME= history -VOLUME_STRUCTURE_FILENAME= direct -VOLUME_ADJ_STRUCTURE_FILENAME= adjoint -SOLUTION_STRUCTURE_FILENAME= direct.dat -RESTART_STRUCTURE_FILENAME= direct.dat -WRT_SOL_FREQ= 9999 -WRT_CON_FREQ= 1 -WRT_CON_FREQ_DUALTIME= 1 -RESTART_SOL= NO +SOLUTION_FILENAME= direct.dat +RESTART_FILENAME= direct.dat +OUTPUT_FILES= PARAVIEW +HISTORY_OUTPUT= COMBO diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 56a34b2f962d..c7e1e907faf3 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -1206,7 +1206,7 @@ def main(): knowlesbeam.cfg_dir = "fea_fsi/MixElemsKnowles" knowlesbeam.cfg_file = "config.cfg" knowlesbeam.test_iter = 0 - knowlesbeam.test_vals = [-14.51360, -13.57735, -28.12642, 9.7306] #last 4 columns + knowlesbeam.test_vals = [-14.51360, -13.57735, -28.12642, 0.44503, 9.7306] #last 5 columns knowlesbeam.su2_exec = "SU2_CFD" knowlesbeam.timeout = 1600 knowlesbeam.tol = 0.0001