From d8ce6e591b4955343c84808af9ec97bc6bf7e5c4 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 10 Apr 2020 23:44:34 +0100 Subject: [PATCH 01/20] CFEASolver::BC_Sym_Plane --- Common/include/CConfig.hpp | 2 +- Common/include/linear_algebra/CSysMatrix.hpp | 8 +- Common/src/linear_algebra/CSysMatrix.cpp | 42 +++++++ SU2_CFD/include/solvers/CFEASolver.hpp | 70 ++++++------ SU2_CFD/include/solvers/CMeshSolver.hpp | 2 +- SU2_CFD/include/solvers/CSolver.hpp | 66 +++++------ .../integration/CStructuralIntegration.cpp | 5 +- SU2_CFD/src/solvers/CFEASolver.cpp | 105 ++++++++++++++---- 8 files changed, 213 insertions(+), 87 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index ab6c1a2f91f6..67030adc3ff5 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -8627,7 +8627,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/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 014355c3952c..97005e4c1eef 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. diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 018c3c84a870..84f62318f831 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() { @@ -1403,6 +1442,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::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double, CSysVector&); template void CSysMatrix::MatrixMatrixAddition(su2double, const CSysMatrix&); #ifdef CODI_REVERSE_TYPE @@ -1413,6 +1453,8 @@ 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::EnforceSolutionAtDOF(unsigned long, unsigned long, passivedouble, CSysVector&); +template void CSysMatrix::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double, CSysVector&); template void CSysMatrix::MatrixMatrixAddition(passivedouble, const CSysMatrix&); template void CSysMatrix::MatrixMatrixAddition(su2double, const CSysMatrix&); #endif diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 02b6f463ae4f..888cc8385e73 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -273,7 +273,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,61 +282,58 @@ 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, @@ -346,48 +343,59 @@ class CFEASolver : public CSolver { /*! * \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_Sym_Plane(CGeometry *geometry, + CNumerics *numerics, + const CConfig *config, + unsigned short val_marker) final; + /*! + * \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. + */ void BC_DispDir(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final; /*! @@ -399,7 +407,7 @@ class CFEASolver : public CSolver { */ inline void BC_Normal_Displacement(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final { } /*! @@ -411,7 +419,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 +431,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 +443,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 +455,7 @@ class CFEASolver : public CSolver { */ void BC_Damper(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final; /*! @@ -459,7 +467,7 @@ class CFEASolver : public CSolver { */ void BC_Deforming(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) final; /*! @@ -841,7 +849,7 @@ class CFEASolver : public CSolver { */ su2double Compute_LoadCoefficient(su2double CurrentTime, su2double RampTime, - CConfig *config) final; + const CConfig *config) final; /*! * \brief A virtual member. 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..c247fba86bb5 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -860,23 +860,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 +884,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 +908,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 +922,7 @@ class CSolver { */ inline virtual void BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! @@ -925,10 +932,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 +947,7 @@ class CSolver { inline virtual void BC_Sine_Load(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! @@ -953,7 +959,7 @@ class CSolver { */ inline virtual void BC_Damper(CGeometry *geometry, CNumerics *numerics, - CConfig *config, + const CConfig *config, unsigned short val_marker) { } /*! @@ -963,10 +969,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) { } /*! @@ -3951,62 +3956,62 @@ 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) { } + 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. */ @@ -4017,13 +4022,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] 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 diff --git a/SU2_CFD/src/integration/CStructuralIntegration.cpp b/SU2_CFD/src/integration/CStructuralIntegration.cpp index ce9b9bd49de8..79eee48be640 100644 --- a/SU2_CFD/src/integration/CStructuralIntegration.cpp +++ b/SU2_CFD/src/integration/CStructuralIntegration.cpp @@ -154,6 +154,9 @@ 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; @@ -163,7 +166,7 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * } } - /*--- Solver linearized system ---*/ + /*--- Solver linear system ---*/ solver_container[MainSolver]->Solve_System(geometry, config); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 5a24ace56ff1..dd1bde87e820 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -849,7 +849,7 @@ void CFEASolver::ResetInitialCondition(CGeometry **geometry, CSolver ***solver_c } -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(); @@ -938,7 +938,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(); @@ -1078,7 +1078,7 @@ 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(); @@ -1158,7 +1158,7 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, C } -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(); @@ -1239,7 +1239,7 @@ void CFEASolver::Compute_MassRes(CGeometry *geometry, CNumerics **numerics, CCon } -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(); @@ -1611,7 +1611,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, } -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. ---*/ @@ -1737,7 +1737,7 @@ 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}; @@ -1769,7 +1769,7 @@ 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(); @@ -1791,7 +1791,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[4] = {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[3] = {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 < 3; ++iDim) + axis = (fabs(normal[iDim]) > fabs(normal[axis]))? iDim : axis; + + if (fabs(normal[axis]) < 0.99*Norm(3,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; @@ -1909,7 +1972,7 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, } -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. ---*/ @@ -2007,7 +2070,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,15 +2078,15 @@ 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 = pow(Load_Dir_Local[0],2) + pow(Load_Dir_Local[1],2); + if (nDim==3) LoadNorm += pow(Load_Dir_Local[2],2); + LoadNorm = sqrt(LoadNorm); 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++) { @@ -2055,7 +2118,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(3,normal); /*--- Compute load vector and update surface load for each node of the boundary element. ---*/ @@ -2070,7 +2133,7 @@ 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)); @@ -2105,7 +2168,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(3,normal); /*--- Compute damping forces. ---*/ @@ -2127,7 +2190,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++) { @@ -2191,7 +2254,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(3,normal) / nNode; /*--- Update area of nodes. ---*/ for (auto iNode = 0u; iNode < nNode; ++iNode) { @@ -2217,7 +2280,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; From c2ad000d8caf062024ce1e4d41993243cdd67f42 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 11 Apr 2020 15:58:07 +0100 Subject: [PATCH 02/20] reduce discrete adjoint memory usage in CMeshSolver --- SU2_CFD/src/solvers/CMeshSolver.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 8fcd4392e01c..9fcc5233713e 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -461,10 +461,16 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig if (multizone) nodes->Set_BGSSolution_k(); - /*--- Compute the stiffness matrix. ---*/ + /*--- 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); - /*--- Clean residual, we do not want an incremental solution. ---*/ + if (ActiveTape) AD::StartRecording(); + + /*--- Clear residual, we do not want an incremental solution. ---*/ SU2_OMP_PARALLEL { LinSysRes.SetValZero(); @@ -499,7 +505,9 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig /*--- 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. ---*/ + 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); From c0e0017fdc0711da3b49f1fe02453add3a748ea5 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sat, 11 Apr 2020 15:58:39 +0100 Subject: [PATCH 03/20] cleanup legacy recording modes --- Common/include/option_structure.hpp | 11 ++---- .../src/drivers/CDiscAdjMultizoneDriver.cpp | 24 ++++++------ .../src/drivers/CDiscAdjSinglezoneDriver.cpp | 10 ++--- SU2_CFD/src/iteration_structure.cpp | 38 +++---------------- 4 files changed, 26 insertions(+), 57 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 4f53145cb13a..41187772993f 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, + CONS_VARS = 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, + COMBINED = 4 }; /*! diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index b5be9b5a904d..a7fe4b95ae33 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(CONS_VARS, Kind_Tape::OBJECTIVE_FUNCTION_TAPE, ZONE_0); RecordingState = NONE; AD::ClearAdjoints(); @@ -212,8 +212,8 @@ 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. + * (2) CONS_VARS: State variables of all solvers in a zone as input. + * (3) MESH_COORDS / MESH_DEFORM: Mesh coordinates as input. * (4) COMBINED: Mesh coordinates and state variables as input. * * By default, all (state and mesh coordinate variables) will be declared as output, @@ -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 != CONS_VARS)) { SetRecording(NONE, Kind_Tape::FULL_TAPE, ZONE_0); - SetRecording(FLOW_CONS_VARS, Kind_Tape::FULL_TAPE, ZONE_0); + SetRecording(CONS_VARS, 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(CONS_VARS, 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 CONS_VARS: 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 == CONS_VARS) { auto solvers = solver_container[iZone][INST_0][MESH_0]; @@ -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 == CONS_VARS)) 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 == CONS_VARS) { 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..9b0312241f9b 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 = CONS_VARS; 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 = CONS_VARS; 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 = CONS_VARS; 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 = CONS_VARS; SecondaryVariables = MESH_COORDS; break; diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index 972cd245e213..4f4e9409fb0c 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -2194,7 +2194,7 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver *****solver, CGeometry ****ge 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 == CONS_VARS || kind_recording == COMBINED) { /*--- Register flow and turbulent variables as input ---*/ @@ -2218,29 +2218,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 ---*/ @@ -2299,8 +2278,7 @@ 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 == COMBINED)){ /*--- Update geometry to get the influence on other geometry variables (normals, volume etc) ---*/ @@ -2670,7 +2648,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, COMBINED); /*--- Run one iteration while tape is passive - this clears all indices ---*/ @@ -2751,9 +2729,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) ---*/ @@ -3270,7 +3245,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 == CONS_VARS || kind_recording == COMBINED){ /*--- Register flow and turbulent variables as input ---*/ @@ -3295,8 +3270,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 == COMBINED)){ /*--- Update geometry to get the influence on other geometry variables (normals, volume etc) ---*/ From 179b5acc570212fc968029aeaf43b70f909bdb18 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 12 Apr 2020 10:54:55 +0100 Subject: [PATCH 04/20] enable DEFORM_LIMIT for CMeshSolver --- SU2_CFD/src/solvers/CMeshSolver.cpp | 41 +++++++++++++++++++++-------- 1 file changed, 30 insertions(+), 11 deletions(-) diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 9fcc5233713e..c2f4f3e81051 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 ---*/ @@ -375,29 +373,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); - + } } } @@ -504,7 +507,7 @@ 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(); @@ -676,6 +679,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){ From c198ae49f0b65c50759f4e3dab0ba502518ac636 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 12 Apr 2020 11:06:03 +0100 Subject: [PATCH 05/20] cleanup unused variables in CSolver and CFEASolver --- SU2_CFD/include/solvers/CFEASolver.hpp | 22 +++++++++------------- SU2_CFD/include/solvers/CSolver.hpp | 2 -- 2 files changed, 9 insertions(+), 15 deletions(-) diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 888cc8385e73..8ad92d9f7b71 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -75,6 +75,15 @@ 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 */ + + 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 */ + #ifdef HAVE_OMP vector > ElemColoring; /*!< \brief Element colors. */ bool LockStrategy = false; /*!< \brief Whether to use an OpenMP lock to guard updates of the Jacobian. */ @@ -182,19 +191,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. */ diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index c247fba86bb5..013189132403 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -173,8 +173,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 */ From 2e9bb598733b462b7fc107b6ffab82443a837a6c Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 12 Apr 2020 15:44:14 +0100 Subject: [PATCH 06/20] compute time residual via numerics for differentiability --- Common/include/linear_algebra/CSysMatrix.hpp | 3 +- Common/src/linear_algebra/CSysMatrix.cpp | 8 +- SU2_CFD/include/solvers/CFEASolver.hpp | 112 +++------ SU2_CFD/include/solvers/CSolver.hpp | 62 ++--- SU2_CFD/src/integration/CIntegration.cpp | 8 +- .../integration/CStructuralIntegration.cpp | 15 +- SU2_CFD/src/iteration_structure.cpp | 12 +- SU2_CFD/src/solvers/CFEASolver.cpp | 217 +++++++++--------- 8 files changed, 183 insertions(+), 254 deletions(-) diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp index 97005e4c1eef..83f6a66ed98d 100644 --- a/Common/include/linear_algebra/CSysMatrix.hpp +++ b/Common/include/linear_algebra/CSysMatrix.hpp @@ -724,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/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp index 84f62318f831..abc2ea6e69b3 100644 --- a/Common/src/linear_algebra/CSysMatrix.cpp +++ b/Common/src/linear_algebra/CSysMatrix.cpp @@ -1368,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. ---*/ @@ -1383,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]; } @@ -1443,7 +1442,6 @@ template void CSysMatrix::InitiateComms(const CSysVector& template void CSysMatrix::CompleteComms(CSysVector&, CGeometry*, CConfig*, unsigned short) const; template void CSysMatrix::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector&); template void CSysMatrix::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double, CSysVector&); -template void CSysMatrix::MatrixMatrixAddition(su2double, const CSysMatrix&); #ifdef CODI_REVERSE_TYPE template class CSysMatrix; @@ -1455,6 +1453,4 @@ template void CSysMatrix::EnforceSolutionAtNode(unsigned long, c template void CSysMatrix::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector&); template void CSysMatrix::EnforceSolutionAtDOF(unsigned long, unsigned long, passivedouble, CSysVector&); template void CSysMatrix::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double, CSysVector&); -template void CSysMatrix::MatrixMatrixAddition(passivedouble, const CSysMatrix&); -template void CSysMatrix::MatrixMatrixAddition(su2double, const CSysMatrix&); #endif diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 8ad92d9f7b71..b05bdf910982 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -79,7 +79,11 @@ class CFEASolver : public CSolver { 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 */ - CSysMatrix MassMatrix; /*!< \brief Sparse structure for storing the mass matrix. */ +#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 */ @@ -172,7 +176,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. @@ -183,7 +187,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, @@ -334,7 +338,7 @@ class CFEASolver : public CSolver { */ void Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, - CConfig *config) final; + const CConfig *config) final; /*! * \brief Compute the dead loads. @@ -394,18 +398,6 @@ class CFEASolver : public CSolver { const CConfig *config, unsigned short val_marker) final; - /*! - * \brief Impose a displacement (constraint) boundary condition. - * \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, - const CConfig *config, - unsigned short val_marker) final { } - /*! * \brief Impose a load boundary condition normal to the boundary. * \param[in] geometry - Geometrical definition of the problem. @@ -473,85 +465,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. @@ -676,45 +643,34 @@ 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, const 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. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. */ - void SetAitken_Relaxation(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) final; + void SetAitken_Relaxation(CGeometry **geometry, const CConfig *config) 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 Update_StructSolution(CGeometry **geometry, CConfig *config) final; /*! * \brief Compute the objective function for a reference geometry @@ -869,6 +825,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/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 013189132403..dd29a0f00758 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -1525,12 +1525,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. @@ -1539,7 +1539,6 @@ class CSolver { * \param[in] config - Definition of the particular problem. */ inline virtual void ImplicitNewmark_Update(CGeometry *geometry, - CSolver **solver_container, CConfig *config) { } /*! @@ -1549,18 +1548,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. @@ -1569,7 +1567,6 @@ class CSolver { * \param[in] config - Definition of the particular problem. */ inline virtual void GeneralizedAlpha_UpdateDisp(CGeometry *geometry, - CSolver **solver_container, CConfig *config) { } /*! @@ -1579,7 +1576,6 @@ class CSolver { * \param[in] config - Definition of the particular problem. */ inline virtual void GeneralizedAlpha_UpdateSolution(CGeometry *geometry, - CSolver **solver_container, CConfig *config) { } /*! @@ -1589,8 +1585,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. @@ -3594,45 +3589,38 @@ class CSolver { /*! * \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. */ - inline virtual void PredictStruct_Displacement(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) { } + inline virtual void PredictStruct_Displacement(CGeometry **geometry, + const 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. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. */ - inline virtual void SetAitken_Relaxation(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) { } + inline virtual void SetAitken_Relaxation(CGeometry **geometry, + const 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 particular problem. */ - inline virtual void Update_StructSolution(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) { } + inline virtual void Update_StructSolution(CGeometry **geometry, + CConfig *config) { } /*! * \brief A virtual member. @@ -4015,7 +4003,7 @@ class CSolver { inline virtual void Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, - CConfig *config) { } + const CConfig *config) { } /*! * \brief A virtual member. diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index d37de0475dab..7f842359b7bd 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -331,12 +331,12 @@ 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); + //if (fsi) solver_container[FEA_SOL]->Update_StructSolution(geometry, 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 79eee48be640..0d0dafe710ea 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; } @@ -160,9 +160,6 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver * 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; } } @@ -174,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; } diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index 4f4e9409fb0c..97d44b5b4476 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -1473,8 +1473,7 @@ void CFEAIteration::Update(COutput *output, /*--- 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]); } } @@ -1496,7 +1495,7 @@ void CFEAIteration::Predictor(COutput *output, /*--- Predict displacements ---*/ - feaSolver->PredictStruct_Displacement(geometry[val_iZone][val_iInst], config[val_iZone], solver[val_iZone][val_iInst]); + feaSolver->PredictStruct_Displacement(geometry[val_iZone][val_iInst], config[val_iZone]); /*--- For parallel simulations we need to communicate the predicted solution before updating the fluid mesh ---*/ @@ -1522,12 +1521,11 @@ 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][INST_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]); + feaSolver->SetAitken_Relaxation(geometry[val_iZone][INST_0], config[val_iZone]); /*----------------- Communicate the predicted solution and the old one ------------------*/ @@ -2825,7 +2823,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(nullptr, config[iZone]); } /*--- MPI dependencies. ---*/ diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index dd1bde87e820..b836e5837820 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1083,6 +1083,10 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, c 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 @@ -1141,10 +1145,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,6 +1160,8 @@ 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, const CConfig *config) { @@ -1163,79 +1169,73 @@ void CFEASolver::Compute_MassRes(CGeometry *geometry, CNumerics **numerics, cons const bool topology_mode = config->GetTopology_Optimization(); const su2double simp_minstiff = config->GetSIMP_MinStiffness(); - /*--- Start OpenMP parallel region. ---*/ + /*--- Clear vector before calculation. ---*/ + TimeRes.SetValZero(); + SU2_OMP_BARRIER - SU2_OMP_PARALLEL - { - /*--- 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]; - 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 } @@ -1330,7 +1330,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(); @@ -1558,7 +1562,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,6 +1613,8 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, } + if (ActiveTape) AD::StartRecording(); + } void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { @@ -1686,7 +1692,7 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, con } -void CFEASolver::Compute_IntegrationConstants(CConfig *config) { +void CFEASolver::Compute_IntegrationConstants(const CConfig *config) { su2double Delta_t= config->GetDelta_DynTime(); @@ -2353,7 +2359,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()); @@ -2413,7 +2419,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'' ---*/ @@ -2429,8 +2435,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; } @@ -2438,7 +2443,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()); @@ -2492,7 +2497,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()); @@ -2556,7 +2561,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()); @@ -2612,7 +2617,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'' ---*/ @@ -2628,8 +2633,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 @@ -2665,7 +2669,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. ---*/ @@ -2681,7 +2685,7 @@ void CFEASolver::GeneralizedAlpha_UpdateDisp(CGeometry *geometry, CSolver **solv } -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); @@ -2740,7 +2744,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(); @@ -2779,14 +2783,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, const 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; @@ -2799,9 +2799,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]; @@ -2809,10 +2809,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]); @@ -2820,7 +2820,7 @@ void CFEASolver::PredictStruct_Displacement(CGeometry **fea_geometry, } break; default: { - fea_nodes->SetSolution_Pred(iPoint); + nodes->SetSolution_Pred(iPoint); } break; } @@ -2828,8 +2828,7 @@ void CFEASolver::PredictStruct_Displacement(CGeometry **fea_geometry, } -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; @@ -2843,7 +2842,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry **fea_geometry, CConfig *fe su2double delta_deltaU[3] = {0.0, 0.0, 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 ---*/ @@ -2854,7 +2853,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) { @@ -2862,8 +2861,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); @@ -2875,10 +2874,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++) { @@ -2920,9 +2919,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, const CConfig *config) { const su2double WAitken = GetWAitken_Dyn(); @@ -2931,14 +2928,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++) { @@ -2948,26 +2945,24 @@ void CFEASolver::SetAitken_Relaxation(CGeometry **fea_geometry, } -void CFEASolver::Update_StructSolution(CGeometry **fea_geometry, - CConfig *fea_config, - CSolver ***fea_solution) { +void CFEASolver::Update_StructSolution(CGeometry **geometry, CConfig *config) { 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); + auto valSolutionPred = nodes->GetSolution_Pred(iPoint); - fea_solution[MESH_0][FEA_SOL]->GetNodes()->SetSolution(iPoint, valSolutionPred); + nodes->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[MESH_0], config, SOLUTION_DISPONLY); + CompleteComms(geometry[MESH_0], config, SOLUTION_DISPONLY); } -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; @@ -3523,7 +3518,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. ---*/ From c0c2b35d36551800de9d71302c791e9fdbed143e Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 12 Apr 2020 18:07:51 +0100 Subject: [PATCH 07/20] remove reduntant MPI comms from CFEASolver --- SU2_CFD/include/solvers/CFEASolver.hpp | 4 +-- SU2_CFD/include/solvers/CSolver.hpp | 13 +------- .../integration/CStructuralIntegration.cpp | 5 ---- SU2_CFD/src/iteration_structure.cpp | 30 +++++-------------- SU2_CFD/src/solvers/CFEASolver.cpp | 12 ++++---- 5 files changed, 16 insertions(+), 48 deletions(-) diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index b05bdf910982..ba768af4cd14 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -338,7 +338,7 @@ class CFEASolver : public CSolver { */ void Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, - const CConfig *config) final; + const CConfig *config); /*! * \brief Compute the dead loads. @@ -646,7 +646,7 @@ class CFEASolver : public CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Configuration of the problem. */ - void PredictStruct_Displacement(CGeometry **geometry, const CConfig *config) final; + void PredictStruct_Displacement(CGeometry **geometry, CConfig *config) final; /*! * \brief Computation of Aitken's coefficient. diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index dd29a0f00758..0352781bb715 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3593,7 +3593,7 @@ class CSolver { * \param[in] config - Definition of the problem. */ inline virtual void PredictStruct_Displacement(CGeometry **geometry, - const CConfig *config) { } + CConfig *config) { } /*! * \brief A virtual member. @@ -3994,17 +3994,6 @@ class CSolver { CNumerics **numerics, const CConfig *config) { } - /*! - * \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. - */ - - inline virtual void Compute_NodalStress(CGeometry *geometry, - CNumerics **numerics, - const CConfig *config) { } - /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/src/integration/CStructuralIntegration.cpp b/SU2_CFD/src/integration/CStructuralIntegration.cpp index 0d0dafe710ea..b73577dba315 100644 --- a/SU2_CFD/src/integration/CStructuralIntegration.cpp +++ b/SU2_CFD/src/integration/CStructuralIntegration.cpp @@ -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 97d44b5b4476..b3ac6de17dce 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -1442,18 +1442,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,9 +1458,9 @@ 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()) + su2double Physical_dt = config[val_iZone]->GetDelta_DynTime(); + 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) { @@ -1493,16 +1487,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]); - /*--- 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); - } + void CFEAIteration::Relaxation(COutput *output, CIntegration ****integration, CGeometry ****geometry, @@ -1553,10 +1541,6 @@ 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], config[val_iZone], config[val_iZone]->GetTimeIter(), diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index b836e5837820..9950b242f0d8 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1910,6 +1910,8 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, bool nonlinear_analysis = (config->GetGeometricConditions() == LARGE_DEFORMATIONS); bool disc_adj_fem = (config->GetKind_Solver() == DISC_ADJ_FEM); + Compute_NodalStress(geometry, numerics, config); + if (disc_adj_fem && nonlinear_analysis) { /*--- For nonlinear discrete adjoint, we have 3 convergence criteria ---*/ @@ -1971,11 +1973,6 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, } - /*--- MPI solution (common to every mode). ---*/ - - InitiateComms(geometry, config, SOLUTION_FEA); - CompleteComms(geometry, config, SOLUTION_FEA); - } void CFEASolver::BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker) { @@ -2783,7 +2780,7 @@ void CFEASolver::Solve_System(CGeometry *geometry, CConfig *config) { } -void CFEASolver::PredictStruct_Displacement(CGeometry **geometry, const CConfig *config) { +void CFEASolver::PredictStruct_Displacement(CGeometry **geometry, CConfig *config) { const unsigned short predOrder = config->GetPredictorOrder(); const su2double Delta_t = config->GetDelta_DynTime(); @@ -2826,6 +2823,9 @@ void CFEASolver::PredictStruct_Displacement(CGeometry **geometry, const CConfig } + InitiateComms(geometry[MESH_0], config, SOLUTION_PRED); + CompleteComms(geometry[MESH_0], config, SOLUTION_PRED); + } void CFEASolver::ComputeAitken_Coefficient(CGeometry **geometry, CConfig *config, unsigned long iOuterIter) { From 2001ea1a2fb85db207ca214770bf566c8ed65c4b Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 12 Apr 2020 18:52:12 +0100 Subject: [PATCH 08/20] more MPI cleanup CMesh/CFEASolver --- Common/include/option_structure.hpp | 1 - SU2_CFD/src/solvers/CFEASolver.cpp | 8 +++--- SU2_CFD/src/solvers/CMeshSolver.cpp | 40 +++++++++++------------------ SU2_CFD/src/solvers/CSolver.cpp | 12 --------- 4 files changed, 19 insertions(+), 42 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 41187772993f..0b245cf188c1 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2082,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/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 9950b242f0d8..059b4fabb49d 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -2677,8 +2677,8 @@ void CFEASolver::GeneralizedAlpha_UpdateDisp(CGeometry *geometry, CConfig *confi /*--- 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); } @@ -2957,8 +2957,8 @@ void CFEASolver::Update_StructSolution(CGeometry **geometry, CConfig *config) { /*--- Perform the MPI communication of the solution, displacements only ---*/ - InitiateComms(geometry[MESH_0], config, SOLUTION_DISPONLY); - CompleteComms(geometry[MESH_0], config, SOLUTION_DISPONLY); + InitiateComms(geometry[MESH_0], config, SOLUTION); + CompleteComms(geometry[MESH_0], config, SOLUTION); } diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index c2f4f3e81051..517f3dcf1e40 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -343,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 { @@ -464,6 +463,16 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig if (multizone) nodes->Set_BGSSolution_k(); + /*--- 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); + + 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(); @@ -473,26 +482,12 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig if (ActiveTape) AD::StartRecording(); - /*--- Clear residual, we do not want an incremental solution. ---*/ + /*--- Clear residual (loses AD info), 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. ---*/ - 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); - /*--- Impose boundary conditions (all of them are ESSENTIAL BC's - displacements). ---*/ SetBoundaryDisplacements(geometry[MESH_0], numerics[FEA_TERM], config); @@ -526,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 ---*/ @@ -540,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); 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]); From 393bae04007dc948ae09b1d9df5327a685846dd4 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 12 Apr 2020 19:10:43 +0100 Subject: [PATCH 09/20] use pre accumulation when computing FE mass matrices --- Common/include/geometry/elements/CElement.hpp | 12 ++++++++++-- SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp | 11 +++++++++-- 2 files changed, 19 insertions(+), 4 deletions(-) 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/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(); + } From 0ed63da6b1a81539f1e513975e510586961a41bc Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 12 Apr 2020 19:23:44 +0100 Subject: [PATCH 10/20] all comms moved to within CFEASolver (less error prone) --- SU2_CFD/include/solvers/CFEASolver.hpp | 2 +- SU2_CFD/include/solvers/CSolver.hpp | 2 +- SU2_CFD/src/iteration_structure.cpp | 5 ----- SU2_CFD/src/solvers/CFEASolver.cpp | 5 ++++- 4 files changed, 6 insertions(+), 8 deletions(-) diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index ba768af4cd14..1af5bf86c2e3 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -663,7 +663,7 @@ class CFEASolver : public CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void SetAitken_Relaxation(CGeometry **geometry, const CConfig *config) final; + void SetAitken_Relaxation(CGeometry **geometry, CConfig *config) final; /*! * \brief Aitken's relaxation of the solution. diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 0352781bb715..7d083f3cb400 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3612,7 +3612,7 @@ class CSolver { * \param[in] config - Definition of the particular problem. */ inline virtual void SetAitken_Relaxation(CGeometry **geometry, - const CConfig *config) { } + CConfig *config) { } /*! * \brief A virtual member. diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index b3ac6de17dce..ec8892ccefe2 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -1515,11 +1515,6 @@ void CFEAIteration::Relaxation(COutput *output, feaSolver->SetAitken_Relaxation(geometry[val_iZone][INST_0], config[val_iZone]); - /*----------------- 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); - } bool CFEAIteration::Monitor(COutput *output, diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 059b4fabb49d..c3d518cd72db 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -2919,7 +2919,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry **geometry, CConfig *config } -void CFEASolver::SetAitken_Relaxation(CGeometry **geometry, const CConfig *config) { +void CFEASolver::SetAitken_Relaxation(CGeometry **geometry, CConfig *config) { const su2double WAitken = GetWAitken_Dyn(); @@ -2943,6 +2943,9 @@ void CFEASolver::SetAitken_Relaxation(CGeometry **geometry, const CConfig *confi } } + InitiateComms(geometry[MESH_0], config, SOLUTION_PRED_OLD); + CompleteComms(geometry[MESH_0], config, SOLUTION_PRED_OLD); + } void CFEASolver::Update_StructSolution(CGeometry **geometry, CConfig *config) { From b2e9f48516b8c097237643b04f713af173e0f7c7 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Sun, 12 Apr 2020 20:11:37 +0100 Subject: [PATCH 11/20] make structural of's accessible for output --- SU2_CFD/include/solvers/CFEASolver.hpp | 27 ++--- SU2_CFD/include/solvers/CSolver.hpp | 16 +-- .../src/drivers/CDiscAdjMultizoneDriver.cpp | 8 +- SU2_CFD/src/iteration_structure.cpp | 14 +-- SU2_CFD/src/output/CElasticityOutput.cpp | 6 +- SU2_CFD/src/solvers/CFEASolver.cpp | 104 ++++-------------- 6 files changed, 49 insertions(+), 126 deletions(-) diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 1af5bf86c2e3..54140c8bed86 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -65,6 +65,7 @@ class CFEASolver : public CSolver { 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. */ @@ -578,6 +579,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 @@ -675,42 +682,30 @@ class CFEASolver : public CSolver { /*! * \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 diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 7d083f3cb400..e1ffbbc9357f 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3687,42 +3687,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) { } + const CConfig *config) { } /*! * \brief A virtual member. diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index a7fe4b95ae33..fa6269aefbcd 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -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; diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index ec8892ccefe2..448ce656d983 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -1398,7 +1398,7 @@ void CFEAIteration::Iterate(COutput *output, } - /*--- Finally, we need to compute the objective function, in case that we are running a discrete adjoint solver... ---*/ + /*--- Finally, we need to compute the objective function, in case we are running the discrete adjoint solver. ---*/ switch (config[val_iZone]->GetKind_ObjFunc()) { case REFERENCE_GEOMETRY: @@ -1406,25 +1406,21 @@ void CFEAIteration::Iterate(COutput *output, 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]); + feaSolver->Compute_OFRefGeom(geometry[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]); + feaSolver->Compute_OFRefNode(geometry[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]); + feaSolver->Compute_OFVolFrac(geometry[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]); + feaSolver->Compute_OFCompliance(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); break; } diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index 0394a6e1eb66..a14f9c92beb6 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -134,6 +134,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 +161,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 +241,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 c3d518cd72db..b0a22ecc6306 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -2081,9 +2081,7 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, const CCo const su2double* Load_Dir_Local = config->GetLoad_Dir(TagBound); /*--- Compute the norm of the vector that was passed in the config file. ---*/ - su2double LoadNorm = pow(Load_Dir_Local[0],2) + pow(Load_Dir_Local[1],2); - if (nDim==3) LoadNorm += pow(Load_Dir_Local[2],2); - LoadNorm = sqrt(LoadNorm); + su2double LoadNorm = Norm(nDim, Load_Dir_Local); su2double CurrentTime=config->GetCurrent_DynTime(); su2double Ramp_Time = config->GetRamp_Time(); @@ -3035,10 +3033,9 @@ void CFEASolver::OutputForwardModeGradient(const 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; @@ -3070,17 +3067,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; @@ -3090,32 +3088,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[3] = {0.0}, dist_reduce[3]; /*--- Convert global point index to local. ---*/ long iPoint = geometry->GetGlobal_to_Local_Point(config->GetRefNode_ID()); @@ -3129,24 +3110,21 @@ void CFEASolver::Compute_OFRefNode(CGeometry *geometry, CSolver **solver_contain SU2_MPI::Allreduce(dist, dist_reduce, 3, 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(3,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; @@ -3156,30 +3134,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 ---*/ @@ -3218,24 +3176,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; + /*--- To be accessible from the output. ---*/ + Total_OFCombo = Total_OFVolFrac; - ofstream myfile_res; - if (config->GetKind_ObjFunc() == TOPOL_DISCRETENESS) - myfile_res.open ("of_topdisc.dat"); - else - myfile_res.open ("of_volfrac.dat"); - - 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(); @@ -3285,17 +3231,9 @@ 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){ From 34d0e7f949416280cc37af6222aa60dac0f60458 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 13 Apr 2020 11:22:36 +0100 Subject: [PATCH 12/20] fix debug build issue with nearest neighbor --- .../interface_interpolation/CNearestNeighbor.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) 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; From fa54eac1405c00c17e2790ddb4a4b12767a62666 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 13 Apr 2020 11:23:04 +0100 Subject: [PATCH 13/20] fix segfault --- SU2_CFD/include/solvers/CFEASolver.hpp | 8 ++++---- SU2_CFD/include/solvers/CSolver.hpp | 8 ++++---- SU2_CFD/src/iteration_structure.cpp | 16 +++++++++------- SU2_CFD/src/solvers/CFEASolver.cpp | 20 ++++++++++---------- 4 files changed, 27 insertions(+), 25 deletions(-) diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 54140c8bed86..95981b44621e 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -653,7 +653,7 @@ class CFEASolver : public CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Configuration of the problem. */ - void PredictStruct_Displacement(CGeometry **geometry, CConfig *config) final; + void PredictStruct_Displacement(CGeometry *geometry, CConfig *config) final; /*! * \brief Computation of Aitken's coefficient. @@ -661,7 +661,7 @@ class CFEASolver : public CSolver { * \param[in] config - Definition of the particular problem. * \param[in] iOuterIter - Current outer iteration. */ - void ComputeAitken_Coefficient(CGeometry **geometry, + void ComputeAitken_Coefficient(CGeometry *geometry, CConfig *config, unsigned long iOuterIter) final; @@ -670,14 +670,14 @@ class CFEASolver : public CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void SetAitken_Relaxation(CGeometry **geometry, CConfig *config) final; + void SetAitken_Relaxation(CGeometry *geometry, CConfig *config) final; /*! * \brief Aitken's relaxation of the solution. * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - void Update_StructSolution(CGeometry **geometry, CConfig *config) final; + void Update_StructSolution(CGeometry *geometry, CConfig *config) final; /*! * \brief Compute the objective function for a reference geometry diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index e1ffbbc9357f..726e6a57a066 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3592,7 +3592,7 @@ class CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the problem. */ - inline virtual void PredictStruct_Displacement(CGeometry **geometry, + inline virtual void PredictStruct_Displacement(CGeometry *geometry, CConfig *config) { } /*! @@ -3601,7 +3601,7 @@ class CSolver { * \param[in] config - Definition of the problem. * \param[in] iOuterIter - Current outer iteration. */ - inline virtual void ComputeAitken_Coefficient(CGeometry **geometry, + inline virtual void ComputeAitken_Coefficient(CGeometry *geometry, CConfig *config, unsigned long iOuterIter) { } @@ -3611,7 +3611,7 @@ class CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - inline virtual void SetAitken_Relaxation(CGeometry **geometry, + inline virtual void SetAitken_Relaxation(CGeometry *geometry, CConfig *config) { } /*! @@ -3619,7 +3619,7 @@ class CSolver { * \param[in] geometry - Geometrical definition of the problem. * \param[in] config - Definition of the particular problem. */ - inline virtual void Update_StructSolution(CGeometry **geometry, + inline virtual void Update_StructSolution(CGeometry *geometry, CConfig *config) { } /*! diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index 448ce656d983..b966a9c88e82 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -1483,7 +1483,7 @@ void CFEAIteration::Predictor(COutput *output, CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL]; - feaSolver->PredictStruct_Displacement(geometry[val_iZone][val_iInst], config[val_iZone]); + feaSolver->PredictStruct_Displacement(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); } @@ -1505,11 +1505,12 @@ void CFEAIteration::Relaxation(COutput *output, /*------------------- Compute the coefficient ---------------------*/ - feaSolver->ComputeAitken_Coefficient(geometry[val_iZone][INST_0], config[val_iZone], 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]); + feaSolver->SetAitken_Relaxation(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); } @@ -1533,7 +1534,8 @@ bool CFEAIteration::Monitor(COutput *output, UsedTime = StopTime - StartTime; 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()); } @@ -1569,13 +1571,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); } @@ -2798,7 +2800,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]); + dir_solver->PredictStruct_Displacement(structural_geometry, config[iZone]); } /*--- MPI dependencies. ---*/ diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index b0a22ecc6306..e3d5c3144851 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -2778,7 +2778,7 @@ void CFEASolver::Solve_System(CGeometry *geometry, CConfig *config) { } -void CFEASolver::PredictStruct_Displacement(CGeometry **geometry, CConfig *config) { +void CFEASolver::PredictStruct_Displacement(CGeometry *geometry, CConfig *config) { const unsigned short predOrder = config->GetPredictorOrder(); const su2double Delta_t = config->GetDelta_DynTime(); @@ -2821,12 +2821,12 @@ void CFEASolver::PredictStruct_Displacement(CGeometry **geometry, CConfig *confi } - InitiateComms(geometry[MESH_0], config, SOLUTION_PRED); - CompleteComms(geometry[MESH_0], config, SOLUTION_PRED); + InitiateComms(geometry, config, SOLUTION_PRED); + CompleteComms(geometry, config, SOLUTION_PRED); } -void CFEASolver::ComputeAitken_Coefficient(CGeometry **geometry, CConfig *config, 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; @@ -2917,7 +2917,7 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry **geometry, CConfig *config } -void CFEASolver::SetAitken_Relaxation(CGeometry **geometry, CConfig *config) { +void CFEASolver::SetAitken_Relaxation(CGeometry *geometry, CConfig *config) { const su2double WAitken = GetWAitken_Dyn(); @@ -2941,12 +2941,12 @@ void CFEASolver::SetAitken_Relaxation(CGeometry **geometry, CConfig *config) { } } - InitiateComms(geometry[MESH_0], config, SOLUTION_PRED_OLD); - CompleteComms(geometry[MESH_0], config, SOLUTION_PRED_OLD); + InitiateComms(geometry, config, SOLUTION_PRED_OLD); + CompleteComms(geometry, config, SOLUTION_PRED_OLD); } -void CFEASolver::Update_StructSolution(CGeometry **geometry, CConfig *config) { +void CFEASolver::Update_StructSolution(CGeometry *geometry, CConfig *config) { SU2_OMP_PARALLEL_(for schedule(static,omp_chunk_size)) for (unsigned long iPoint=0; iPoint < nPointDomain; iPoint++) { @@ -2958,8 +2958,8 @@ void CFEASolver::Update_StructSolution(CGeometry **geometry, CConfig *config) { /*--- Perform the MPI communication of the solution, displacements only ---*/ - InitiateComms(geometry[MESH_0], config, SOLUTION); - CompleteComms(geometry[MESH_0], config, SOLUTION); + InitiateComms(geometry, config, SOLUTION); + CompleteComms(geometry, config, SOLUTION); } From d3a693368486b21792ab7310af8dfcd7144617a5 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 13 Apr 2020 12:47:54 +0100 Subject: [PATCH 14/20] compute structural of's every iteration for monitoring --- SU2_CFD/include/solvers/CFEASolver.hpp | 2 +- SU2_CFD/include/solvers/CSolver.hpp | 13 --- SU2_CFD/src/iteration_structure.cpp | 26 ------ SU2_CFD/src/solvers/CFEASolver.cpp | 87 ++++++++++++++------ SU2_CFD/src/solvers/CRadP1Solver.cpp | 1 - TestCases/fea_fsi/MixElemsKnowles/config.cfg | 2 +- TestCases/serial_regression.py | 2 +- 7 files changed, 63 insertions(+), 70 deletions(-) diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 95981b44621e..05b90b52ae76 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -717,7 +717,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. diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 726e6a57a066..004b2a764ea8 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. */ @@ -3716,18 +3715,6 @@ class CSolver { inline virtual void Compute_OFCompliance(CGeometry *geometry, 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] 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) { } - /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index b966a9c88e82..c7129a4a2a15 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -1398,32 +1398,6 @@ void CFEAIteration::Iterate(COutput *output, } - /*--- Finally, we need to compute the objective function, in case we are running the 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], 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], config[val_iZone]); - break; - case VOLUME_FRACTION: - case TOPOL_DISCRETENESS: - feaSolver->Compute_OFVolFrac(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); - break; - case TOPOL_COMPLIANCE: - feaSolver->Compute_OFCompliance(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); - break; - } - } void CFEAIteration::Update(COutput *output, diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index e3d5c3144851..013fe3cb2eee 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -188,9 +188,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. ---*/ @@ -1904,14 +1901,31 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, CNumerics *numerics, const CCon 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); + const bool disc_adj_fem = (config->GetKind_Solver() == DISC_ADJ_FEM); - bool nonlinear_analysis = (config->GetGeometricConditions() == LARGE_DEFORMATIONS); - bool disc_adj_fem = (config->GetKind_Solver() == DISC_ADJ_FEM); + /*--- Compute stresses for monitoring and output. ---*/ Compute_NodalStress(geometry, numerics, config); + /*--- Compute the objective function to be able to monitor it. ---*/ + + const auto kindObjFunc = config->GetKind_ObjFunc(); + + if (((kindObjFunc == REFERENCE_GEOMETRY) || (kindObjFunc == REFERENCE_NODE)) && + ((config->GetDV_FEA() == YOUNG_MODULUS) || (config->GetDV_FEA() == DENSITY_VAL))) { + + Stiffness_Penalty(geometry, solver_container, numerics, config); + } + + 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 (disc_adj_fem && nonlinear_analysis) { /*--- For nonlinear discrete adjoint, we have 3 convergence criteria ---*/ @@ -1931,46 +1945,59 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, /*--- If the problem is linear, the only check we do is the RMS of the residuals. ---*/ /*--- Compute the residual Ax-f ---*/ + CSysVector LinSysAux(nPoint, nPointDomain, nVar, nullptr); + 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); + } // end SU2_OMP_PARALLEL + } } @@ -3238,6 +3265,12 @@ void CFEASolver::Compute_OFCompliance(CGeometry *geometry, const CConfig *config 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; 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/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/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 From 5ffabccdbd030ec7437e838eec865c894c653c10 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 13 Apr 2020 13:18:19 +0100 Subject: [PATCH 15/20] fix directdiff build --- SU2_CFD/src/solvers/CFEASolver.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 013fe3cb2eee..a2fdc45f93d8 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1945,7 +1945,11 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, /*--- 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 { From 760531cb57fd9b6be389a0b3a3e0f1f2caac6680 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 13 Apr 2020 18:47:35 +0100 Subject: [PATCH 16/20] fix topology optimization script and example --- SU2_CFD/src/output/CAdjElasticityOutput.cpp | 18 +++++++++------ SU2_CFD/src/solvers/CFEASolver.cpp | 21 +++++++---------- SU2_PY/topology_optimization.py | 23 +++++++++---------- .../fea_topology/quick_start/settings.cfg | 22 +++++------------- ...gs_refnode.cfg => settings_compliance.cfg} | 23 ++++++------------- .../quick_start/settings_volfrac.cfg | 22 +++++++----------- 6 files changed, 51 insertions(+), 78 deletions(-) rename TestCases/fea_topology/quick_start/{settings_refnode.cfg => settings_compliance.cfg} (66%) 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/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index a2fdc45f93d8..202465c3b3f8 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -2029,13 +2029,8 @@ void CFEASolver::BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, const 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[4] = {0}; + su2double nodeCoord_ref[4][3] = {{0.0}}, nodeCoord_curr[4][3] = {{0.0}}; /*--- Identify the kind of boundary element. ---*/ @@ -2061,8 +2056,8 @@ void CFEASolver::BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, const /*--- 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[3] = {0.0}; + su2double normal_curr[3] = {0.0}; switch (nNodes) { case 2: LineNormal(nodeCoord_ref, normal_ref); break; @@ -2124,9 +2119,9 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, const CCo 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[4] = {0}; - const su2double* nodeCoord[4] = {nullptr, nullptr, nullptr, nullptr}; + const su2double* nodeCoord[4] = {nullptr}; /*--- Identify the kind of boundary element. ---*/ @@ -2142,7 +2137,7 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, const CCo /*--- Compute area of the boundary element. ---*/ - su2double normal[3] = {0.0, 0.0, 0.0}; + su2double normal[3] = {0.0}; switch (nNodes) { case 2: LineNormal(nodeCoord, normal); break; @@ -2230,7 +2225,7 @@ void CFEASolver::BC_Deforming(CGeometry *geometry, CNumerics *numerics, const CC auto iNode = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Retrieve the boundary displacement ---*/ - su2double Disp[3] = {0.0, 0.0, 0.0}; + su2double Disp[3] = {0.0}; for (unsigned short iDim = 0; iDim < nDim; iDim++) Disp[iDim] = nodes->GetBound_Disp(iNode,iDim); 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_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 From 94a6cd438a4929d7d5286f8358eb8a67730b6c4f Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 14 Apr 2020 12:09:02 +0100 Subject: [PATCH 17/20] improve the incremental load approach for multizone cases --- Common/include/CConfig.hpp | 20 ++---- SU2_CFD/include/solvers/CFEASolver.hpp | 31 +++----- SU2_CFD/include/solvers/CSolver.hpp | 19 ++--- SU2_CFD/include/variables/CFEAVariable.hpp | 4 +- SU2_CFD/include/variables/CVariable.hpp | 2 +- SU2_CFD/src/drivers/CMultizoneDriver.cpp | 18 +++-- SU2_CFD/src/iteration_structure.cpp | 83 +++++++++------------- SU2_CFD/src/output/CElasticityOutput.cpp | 7 +- SU2_CFD/src/solvers/CFEASolver.cpp | 60 +++++++++------- 9 files changed, 102 insertions(+), 142 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 67030adc3ff5..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; } /*! diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 05b90b52ae76..48ba1e6e839a 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -48,18 +48,18 @@ 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. */ @@ -242,18 +242,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. @@ -760,7 +748,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 diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 004b2a764ea8..cb1240d2765f 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3574,18 +3574,6 @@ class CSolver { CConfig *config, unsigned long TimeIter) { } - /*! - * \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] geometry - Geometrical definition of the problem. @@ -4025,15 +4013,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/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/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index c7129a4a2a15..f6955c9da628 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; - - Criteria_UTOL = config[val_iZone]->GetIncLoad_Criteria(0); - Criteria_RTOL = config[val_iZone]->GetIncLoad_Criteria(1); - Criteria_ETOL = config[val_iZone]->GetIncLoad_Criteria(2); + /*--- Early return if we already meet the convergence criteria. ---*/ + if (StopCalc) return; - Residual_UTOL = log10(feaSolver->LinSysSol.norm()); - Residual_RTOL = log10(feaSolver->LinSysRes.norm()); - Residual_ETOL = log10(feaSolver->LinSysSol.dot(feaSolver->LinSysRes)); + /*--- Check user-defined criteria to either increment loads or continue with NR iterations. ---*/ - meetCriteria = ( ( Residual_UTOL < Criteria_UTOL ) && - ( Residual_RTOL < Criteria_RTOL ) && - ( Residual_ETOL < Criteria_ETOL ) ); + bool meetCriteria = true; + for (int i = 0; i < 3; ++i) + meetCriteria &= (feaSolver->GetRes_FEM(i) < config[val_iZone]->GetIncLoad_Criteria(i)); - /*--- 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 ---*/ @@ -1393,7 +1377,8 @@ void CFEAIteration::Iterate(COutput *output, } } - + /*--- Just to be sure, set default increment settings. ---*/ + feaSolver->SetLoad_Increment(0, 1.0); } } diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index a14f9c92beb6..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){ diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 202465c3b3f8..eb5ec7df1861 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); @@ -823,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; @@ -832,18 +834,20 @@ 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, const CConfig *config) { @@ -1902,7 +1906,6 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, CNumerics **numerics, unsigned short iMesh) { const bool nonlinear_analysis = (config->GetGeometricConditions() == LARGE_DEFORMATIONS); - const bool disc_adj_fem = (config->GetKind_Solver() == DISC_ADJ_FEM); /*--- Compute stresses for monitoring and output. ---*/ @@ -1926,21 +1929,28 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CSolver **solver_container, case TOPOL_COMPLIANCE: Compute_OFCompliance(geometry, config); break; } - if (disc_adj_fem && nonlinear_analysis) { + if (nonlinear_analysis) { - /*--- For nonlinear discrete adjoint, we have 3 convergence criteria ---*/ + /*--- 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 ---*/ - /*--- 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 -----------*/ - - 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 + 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 } - - if (!nonlinear_analysis) { + else { /*--- If the problem is linear, the only check we do is the RMS of the residuals. ---*/ /*--- Compute the residual Ax-f ---*/ From 14963358bacde10748fa61041e61f1e4cf7c0edb Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 14 Apr 2020 12:25:19 +0100 Subject: [PATCH 18/20] small fix --- SU2_CFD/src/iteration_structure.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index f6955c9da628..0f89df7c7351 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -1310,7 +1310,7 @@ void CFEAIteration::Iterate(COutput *output, bool meetCriteria = true; for (int i = 0; i < 3; ++i) - meetCriteria &= (feaSolver->GetRes_FEM(i) < config[val_iZone]->GetIncLoad_Criteria(i)); + meetCriteria &= (log10(feaSolver->GetRes_FEM(i)) < config[val_iZone]->GetIncLoad_Criteria(i)); /*--- If the criteria is met, i.e. the load is not too large, continue the regular calculation. ---*/ @@ -1368,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); From af600225ecd6ff5bfe88a5e0d5e41f57341efd09 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 14 Apr 2020 19:58:05 +0100 Subject: [PATCH 19/20] address review comments --- Common/include/option_structure.hpp | 4 +- SU2_CFD/include/iteration_structure.hpp | 2 +- SU2_CFD/include/solvers/CFEASolver.hpp | 11 +- SU2_CFD/include/solvers/CSolver.hpp | 9 -- .../src/drivers/CDiscAdjMultizoneDriver.cpp | 20 ++-- .../src/drivers/CDiscAdjSinglezoneDriver.cpp | 8 +- SU2_CFD/src/integration/CIntegration.cpp | 1 - SU2_CFD/src/iteration_structure.cpp | 14 +-- SU2_CFD/src/solvers/CFEASolver.cpp | 105 +++++++----------- 9 files changed, 69 insertions(+), 105 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 0b245cf188c1..725167f7d878 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2020,10 +2020,10 @@ static const MapType DirectDiff_Var_Map = { enum ENUM_RECORDING { - CONS_VARS = 1, + SOLUTION_VARIABLES = 1, MESH_COORDS = 2, MESH_DEFORM = 3, - COMBINED = 4 + SOLUTION_AND_MESH = 4 }; /*! diff --git a/SU2_CFD/include/iteration_structure.hpp b/SU2_CFD/include/iteration_structure.hpp index 5291eb00380f..f099ba880fef 100644 --- a/SU2_CFD/include/iteration_structure.hpp +++ b/SU2_CFD/include/iteration_structure.hpp @@ -1195,7 +1195,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, diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 48ba1e6e839a..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}; @@ -660,13 +662,6 @@ class CFEASolver : public CSolver { */ void SetAitken_Relaxation(CGeometry *geometry, CConfig *config) final; - /*! - * \brief Aitken's relaxation of the solution. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - void Update_StructSolution(CGeometry *geometry, CConfig *config) final; - /*! * \brief Compute the objective function for a reference geometry * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index cb1240d2765f..bc5b71a6ffbb 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -3592,7 +3592,6 @@ class CSolver { CConfig *config, unsigned long iOuterIter) { } - /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. @@ -3601,14 +3600,6 @@ class CSolver { inline virtual void SetAitken_Relaxation(CGeometry *geometry, CConfig *config) { } - /*! - * \brief A virtual member. - * \param[in] geometry - Geometrical definition of the problem. - * \param[in] config - Definition of the particular problem. - */ - inline virtual void Update_StructSolution(CGeometry *geometry, - CConfig *config) { } - /*! * \brief A virtual member. * \param[in] geometry - Geometrical definition of the problem. diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index fa6269aefbcd..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(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) CONS_VARS: State variables of all solvers in a zone 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) COMBINED: Mesh coordinates and state variables 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 != CONS_VARS)) { + if (!retape && (RecordingState != SOLUTION_VARIABLES)) { SetRecording(NONE, Kind_Tape::FULL_TAPE, ZONE_0); - SetRecording(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(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. ---*/ @@ -430,7 +430,7 @@ void CDiscAdjMultizoneDriver::SetRecording(unsigned short kind_recording, Kind_T 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 CONS_VARS: cout << "Storing computational graph wrt CONSERVATIVE VARIABLES." << 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 == CONS_VARS) { + if (rank == MASTER_NODE && kind_recording == SOLUTION_VARIABLES) { auto solvers = solver_container[iZone][INST_0][MESH_0]; @@ -762,14 +762,14 @@ void CDiscAdjMultizoneDriver::SetObjFunction(unsigned short kind_recording) { break; } - if (ObjectiveNotCovered && (rank == MASTER_NODE) && (kind_recording == 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 == 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 9b0312241f9b..3e068ebd4726 100644 --- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp @@ -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 = 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 = 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 = CONS_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 = CONS_VARS; + MainVariables = SOLUTION_VARIABLES; SecondaryVariables = MESH_COORDS; break; diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index 7f842359b7bd..5b72073a9546 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -334,7 +334,6 @@ void CIntegration::SetStructural_Solver(CGeometry *geometry, CSolver **solver_co if (fsi) solver_container[FEA_SOL]->ImplicitNewmark_Relaxation(geometry, config); break; case (GENERALIZED_ALPHA): - //if (fsi) solver_container[FEA_SOL]->Update_StructSolution(geometry, config); solver_container[FEA_SOL]->GeneralizedAlpha_UpdateSolution(geometry, config); solver_container[FEA_SOL]->GeneralizedAlpha_UpdateLoads(geometry, config); break; diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index 0f89df7c7351..40cf7b13ad48 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -1412,8 +1412,8 @@ void CFEAIteration::Update(COutput *output, /*--- Verify convergence criteria (based on total time) ---*/ - su2double Physical_dt = config[val_iZone]->GetDelta_DynTime(); - su2double Physical_t = (TimeIter+1)*Physical_dt; + 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); @@ -2127,7 +2127,7 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver *****solver, CGeometry ****ge bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); bool heat = config[iZone]->GetWeakly_Coupled_Heat(); - if (kind_recording == CONS_VARS || kind_recording == COMBINED) { + if (kind_recording == SOLUTION_VARIABLES || kind_recording == SOLUTION_AND_MESH) { /*--- Register flow and turbulent variables as input ---*/ @@ -2211,7 +2211,7 @@ 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 == COMBINED)){ + 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) ---*/ @@ -2581,7 +2581,7 @@ void CDiscAdjFEAIteration::SetRecording(COutput *output, /*--- Clear indices of coupling variables ---*/ - SetDependencies(solver, geometry, numerics, config, val_iZone, val_iInst, COMBINED); + SetDependencies(solver, geometry, numerics, config, val_iZone, val_iInst, SOLUTION_AND_MESH); /*--- Run one iteration while tape is passive - this clears all indices ---*/ @@ -3178,7 +3178,7 @@ void CDiscAdjHeatIteration::RegisterInput(CSolver *****solver, unsigned short iZone, unsigned short iInst, unsigned short kind_recording){ - if (kind_recording == CONS_VARS || kind_recording == COMBINED){ + if (kind_recording == SOLUTION_VARIABLES || kind_recording == SOLUTION_AND_MESH){ /*--- Register flow and turbulent variables as input ---*/ @@ -3203,7 +3203,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 == COMBINED)){ + 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/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index eb5ec7df1861..1c1864f5f667 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -529,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; @@ -682,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]); @@ -885,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++) { @@ -978,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++) { @@ -1116,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); @@ -1195,7 +1195,7 @@ void CFEASolver::Compute_MassRes(CGeometry *geometry, CNumerics **numerics, cons 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); @@ -1277,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++) { @@ -1382,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++) { @@ -1650,7 +1650,7 @@ void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, con 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); @@ -1747,7 +1747,7 @@ void CFEASolver::Compute_IntegrationConstants(const CConfig *config) { 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++) { @@ -1780,7 +1780,7 @@ void CFEASolver::BC_Clamped_Post(CGeometry *geometry, CNumerics *numerics, const 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++) { @@ -1806,7 +1806,7 @@ void CFEASolver::BC_Sym_Plane(CGeometry *geometry, CNumerics *numerics, const CC /*--- Determine axis of symmetry based on the normal of the first element in the marker. ---*/ - const su2double* nodeCoord[4] = {nullptr}; + const su2double* nodeCoord[MAXNNODE_2D] = {nullptr}; const bool quad = (geometry->bound[val_marker][0]->GetVTK_Type() == QUADRILATERAL); const unsigned short nNodes = quad? 4 : nDim; @@ -1816,7 +1816,7 @@ void CFEASolver::BC_Sym_Plane(CGeometry *geometry, CNumerics *numerics, const CC nodeCoord[iNode] = geometry->node[iPoint]->GetCoord(); } - su2double normal[3] = {0.0}; + su2double normal[MAXNDIM] = {0.0}; switch (nNodes) { case 2: LineNormal(nodeCoord, normal); break; @@ -1825,10 +1825,10 @@ void CFEASolver::BC_Sym_Plane(CGeometry *geometry, CNumerics *numerics, const CC } auto axis = 0u; - for (auto iDim = 1u; iDim < 3; ++iDim) + for (auto iDim = 1u; iDim < MAXNDIM; ++iDim) axis = (fabs(normal[iDim]) > fabs(normal[axis]))? iDim : axis; - if (fabs(normal[axis]) < 0.99*Norm(3,normal)) { + if (fabs(normal[axis]) < 0.99*Norm(int(MAXNDIM),normal)) { SU2_MPI::Error("The structural solver only supports axis-aligned symmetry planes.",CURRENT_FUNCTION); } @@ -1869,13 +1869,7 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, CNumerics *numerics, const CCon 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(); @@ -1883,6 +1877,7 @@ void CFEASolver::BC_DispDir(CGeometry *geometry, CNumerics *numerics, const CCon su2double TotalDisp = ModAmpl * DispDirVal * DispDirMult / DispDirMod; + su2double DispDir[MAXNVAR] = {0.0}; for (iDim = 0; iDim < nDim; iDim++) DispDir[iDim] = TotalDisp * DispDirLocal[iDim]; @@ -2039,8 +2034,9 @@ void CFEASolver::BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, const for (unsigned long iElem = 0; iElem < geometry->GetnElem_Bound(val_marker); iElem++) { unsigned short iNode, iDim; - unsigned long indexNode[4] = {0}; - su2double nodeCoord_ref[4][3] = {{0.0}}, nodeCoord_curr[4][3] = {{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. ---*/ @@ -2066,8 +2062,8 @@ void CFEASolver::BC_Normal_Load(CGeometry *geometry, CNumerics *numerics, const /*--- Compute area vectors in reference and current configurations. ---*/ - su2double normal_ref[3] = {0.0}; - su2double normal_curr[3] = {0.0}; + su2double normal_ref[MAXNDIM] = {0.0}; + su2double normal_curr[MAXNDIM] = {0.0}; switch (nNodes) { case 2: LineNormal(nodeCoord_ref, normal_ref); break; @@ -2129,9 +2125,9 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, const CCo for (unsigned long iElem = 0; iElem < geometry->GetnElem_Bound(val_marker); iElem++) { unsigned short iNode, iDim; - unsigned long indexNode[4] = {0}; + unsigned long indexNode[MAXNNODE_2D] = {0}; - const su2double* nodeCoord[4] = {nullptr}; + const su2double* nodeCoord[MAXNNODE_2D] = {nullptr}; /*--- Identify the kind of boundary element. ---*/ @@ -2147,7 +2143,7 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, const CCo /*--- Compute area of the boundary element. ---*/ - su2double normal[3] = {0.0}; + su2double normal[MAXNDIM] = {0.0}; switch (nNodes) { case 2: LineNormal(nodeCoord, normal); break; @@ -2155,7 +2151,7 @@ void CFEASolver::BC_Dir_Load(CGeometry *geometry, CNumerics *numerics, const CCo case 4: QuadrilateralNormal(nodeCoord, normal); break; } - su2double area = Norm(3,normal); + su2double area = Norm(int(MAXNDIM),normal); /*--- Compute load vector and update surface load for each node of the boundary element. ---*/ @@ -2177,9 +2173,9 @@ void CFEASolver::BC_Damper(CGeometry *geometry, CNumerics *numerics, const CConf 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; @@ -2197,7 +2193,7 @@ void CFEASolver::BC_Damper(CGeometry *geometry, CNumerics *numerics, const CConf /*--- 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; @@ -2205,7 +2201,7 @@ void CFEASolver::BC_Damper(CGeometry *geometry, CNumerics *numerics, const CConf case 4: QuadrilateralNormal(nodeCoord, normal); break; } - su2double area = Norm(3,normal); + su2double area = Norm(int(MAXNDIM),normal); /*--- Compute damping forces. ---*/ @@ -2235,7 +2231,7 @@ void CFEASolver::BC_Deforming(CGeometry *geometry, CNumerics *numerics, const CC auto iNode = geometry->vertex[val_marker][iVertex]->GetNode(); /*--- Retrieve the boundary displacement ---*/ - su2double Disp[3] = {0.0}; + su2double Disp[MAXNVAR] = {0.0}; for (unsigned short iDim = 0; iDim < nDim; iDim++) Disp[iDim] = nodes->GetBound_Disp(iNode,iDim); @@ -2270,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; @@ -2283,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; @@ -2291,7 +2287,7 @@ void CFEASolver::Integrate_FSI_Loads(CGeometry *geometry, const CConfig *config) case 4: QuadrilateralNormal(coords, normal); break; } - su2double area = Norm(3,normal) / nNode; + su2double area = Norm(int(MAXNDIM),normal) / nNode; /*--- Update area of nodes. ---*/ for (auto iNode = 0u; iNode < nNode; ++iNode) { @@ -2309,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); @@ -2872,8 +2868,8 @@ void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, CConfig *config, 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 = config->GetRelaxation_Method_FSI(); @@ -2982,23 +2978,6 @@ void CFEASolver::SetAitken_Relaxation(CGeometry *geometry, CConfig *config) { } -void CFEASolver::Update_StructSolution(CGeometry *geometry, CConfig *config) { - - SU2_OMP_PARALLEL_(for schedule(static,omp_chunk_size)) - for (unsigned long iPoint=0; iPoint < nPointDomain; iPoint++) { - - auto valSolutionPred = nodes->GetSolution_Pred(iPoint); - - nodes->SetSolution(iPoint, valSolutionPred); - } - - /*--- Perform the MPI communication of the solution, displacements only ---*/ - - InitiateComms(geometry, config, SOLUTION); - CompleteComms(geometry, config, SOLUTION); - -} - void CFEASolver::OutputForwardModeGradient(const CConfig *config, bool newFile, su2double fun, su2double fun_avg, su2double der, su2double der_avg) const { @@ -3132,7 +3111,7 @@ void CFEASolver::Compute_OFRefNode(CGeometry *geometry, const CConfig *config){ bool fsi = config->GetFSI_Simulation(); unsigned long TimeIter = config->GetTimeIter(); - su2double dist[3] = {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()); @@ -3144,9 +3123,9 @@ void CFEASolver::Compute_OFRefNode(CGeometry *geometry, const CConfig *config){ } } - 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 = config->GetRefNode_Penalty() * Norm(3,dist_reduce) + PenaltyValue; + Total_OFRefNode = config->GetRefNode_Penalty() * Norm(int(MAXNVAR),dist_reduce) + PenaltyValue; Global_OFRefNode += Total_OFRefNode; @@ -3297,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); From 4b810468db36745aac1276f2fc6433285c26512f Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Thu, 16 Apr 2020 23:19:19 +0100 Subject: [PATCH 20/20] remove one more obsolete method for discadj fsi --- SU2_CFD/include/iteration_structure.hpp | 26 ----- SU2_CFD/src/iteration_structure.cpp | 130 +++++++----------------- 2 files changed, 35 insertions(+), 121 deletions(-) diff --git a/SU2_CFD/include/iteration_structure.hpp b/SU2_CFD/include/iteration_structure.hpp index f099ba880fef..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, @@ -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/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp index 40cf7b13ad48..54efa417df01 100644 --- a/SU2_CFD/src/iteration_structure.cpp +++ b/SU2_CFD/src/iteration_structure.cpp @@ -2048,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]); } @@ -2123,7 +2110,6 @@ 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(); @@ -2131,9 +2117,7 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver *****solver, CGeometry ****ge /*--- 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]); @@ -2176,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]); } } @@ -2211,6 +2197,7 @@ 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 == SOLUTION_AND_MESH)){ /*--- Update geometry to get the influence on other geometry variables (normals, volume etc) ---*/ @@ -2248,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]); @@ -2277,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, @@ -2330,6 +2282,7 @@ void CDiscAdjFluidIteration::Update(COutput *output, } } } + bool CDiscAdjFluidIteration::Monitor(COutput *output, CIntegration ****integration, CGeometry ****geometry, @@ -2804,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,