From ec4eb9fc013d0c0422b70603719c31d34c3c3c9e Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 19 Jan 2021 15:47:02 +0000 Subject: [PATCH 1/9] initial implementation, not working with AD --- Common/include/CConfig.hpp | 8 +- Common/include/geometry/elements/CElement.hpp | 7 ++ Common/include/option_structure.hpp | 4 +- Common/src/CConfig.cpp | 7 +- SU2_CFD/include/solvers/CFEASolver.hpp | 11 +- SU2_CFD/include/solvers/CSolver.hpp | 10 +- .../numerics/elasticity/CFEAElasticity.cpp | 9 ++ .../elasticity/CFEALinearElasticity.cpp | 14 ++- .../elasticity/CFEANonlinearElasticity.cpp | 19 ++- SU2_CFD/src/output/CElasticityOutput.cpp | 3 + SU2_CFD/src/solvers/CFEASolver.cpp | 112 +++++++++++------- 11 files changed, 146 insertions(+), 58 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 9682fd95ad3e..d840895410c5 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -990,7 +990,8 @@ class CConfig { bool FEAAdvancedMode; /*!< \brief Determine if advanced features are used from the element-based FEA analysis (experimental). */ 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. */ + DV_Penalty, /*!< \brief Penalty weight to add a constraint to the total amount of stiffness. */ + AllowedVMStress; /*!< \brief Maximum stress for the stress penalty objective function. */ 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. */ su2double ParMETIS_tolerance; /*!< \brief Load balancing tolerance for ParMETIS. */ @@ -8532,6 +8533,11 @@ class CConfig { */ su2double GetTotalDV_Penalty(void) const { return DV_Penalty; } + /*! + * \brief Get the maximum allowed VM stress for the stress penalty objective function. + */ + su2double GetAllowedVMStress(void) const { return AllowedVMStress; } + /*! * \brief Get whether a predictor is used for FSI applications. * \return Bool: determines if predictor is used or not diff --git a/Common/include/geometry/elements/CElement.hpp b/Common/include/geometry/elements/CElement.hpp index 8e86e8ece7c8..543023fbae98 100644 --- a/Common/include/geometry/elements/CElement.hpp +++ b/Common/include/geometry/elements/CElement.hpp @@ -470,6 +470,13 @@ class CElement { AD::SetPreaccOut(Mab.data(), nNodes*nNodes); } + /*! + * \brief Register the dead load as a pre-accumulation output. + */ + inline void SetPreaccOut_FDL_a(void) { + AD::SetPreaccOut(FDL_a.data(), nNodes*nDim); + } + }; /*! diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index c69c5f5d44ca..e39eba865c34 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1524,7 +1524,8 @@ enum ENUM_OBJECTIVE { REFERENCE_NODE = 61, /*!< \brief Objective function defined as the difference of a particular node respect to a reference position. */ VOLUME_FRACTION = 62, /*!< \brief Volume average physical density, for material-based topology optimization applications. */ TOPOL_DISCRETENESS = 63, /*!< \brief Measure of the discreteness of the current topology. */ - TOPOL_COMPLIANCE = 64 /*!< \brief Measure of the discreteness of the current topology. */ + TOPOL_COMPLIANCE = 64, /*!< \brief Measure of the discreteness of the current topology. */ + STRESS_PENALTY = 65, /*!< \brief Penalty function of VM stresses above a maximum value. */ }; static const MapType Objective_Map = { MakePair("DRAG", DRAG_COEFFICIENT) @@ -1575,6 +1576,7 @@ static const MapType Objective_Map = { MakePair("VOLUME_FRACTION", VOLUME_FRACTION) MakePair("TOPOL_DISCRETENESS", TOPOL_DISCRETENESS) MakePair("TOPOL_COMPLIANCE", TOPOL_COMPLIANCE) + MakePair("STRESS_PENALTY", STRESS_PENALTY) }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index dc370d2e843d..5364c48430f2 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2289,11 +2289,14 @@ void CConfig::SetConfig_Options() { /*!\brief REFERENCE_NODE\n DESCRIPTION: Reference node for the structure (optimization applications) */ addUnsignedLongOption("REFERENCE_NODE", refNodeID, 0); - /* DESCRIPTION: Modulus of the electric fields */ + /*!\brief REFERENCE_NODE_DISPLACEMENT\n DESCRIPTION: Target displacement of the reference node \ingroup Config*/ addDoubleListOption("REFERENCE_NODE_DISPLACEMENT", nDim_RefNode, RefNode_Displacement); /*!\brief REFERENCE_NODE_PENALTY\n DESCRIPTION: Penalty weight value for the objective function \ingroup Config*/ addDoubleOption("REFERENCE_NODE_PENALTY", RefNode_Penalty, 1E3); + /*!\brief ALLOWED_VONMISSES_STRESS\n DESCRIPTION: Maximum allowed stress for structural optimization \ingroup Config*/ + addDoubleOption("ALLOWED_VONMISSES_STRESS", AllowedVMStress, 1.0); + /*!\brief REGIME_TYPE \n DESCRIPTION: Geometric condition \n OPTIONS: see \link Struct_Map \endlink \ingroup Config*/ addEnumOption("GEOMETRIC_CONDITIONS", Kind_Struct_Solver, Struct_Map, SMALL_DEFORMATIONS); /*!\brief REGIME_TYPE \n DESCRIPTION: Material model \n OPTIONS: see \link Material_Map \endlink \ingroup Config*/ @@ -6051,6 +6054,7 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) { case VOLUME_FRACTION: cout << "Volume fraction objective function." << endl; break; case TOPOL_DISCRETENESS: cout << "Topology discreteness objective function." << endl; break; case TOPOL_COMPLIANCE: cout << "Topology compliance objective function." << endl; break; + case STRESS_PENALTY: cout << "Stress penalty objective function." << endl; break; } } else { @@ -8018,6 +8022,7 @@ string CConfig::GetObjFunc_Extension(string val_filename) const { case VOLUME_FRACTION: AdjExt = "_volfrac"; break; case TOPOL_DISCRETENESS: AdjExt = "_topdisc"; break; case TOPOL_COMPLIANCE: AdjExt = "_topcomp"; break; + case STRESS_PENALTY: AdjExt = "_stress"; break; } } else{ diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 31ec0790bad5..757eafae6798 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -69,6 +69,7 @@ class CFEASolver : public CSolver { su2double Total_OFVolFrac; /*!< \brief Total Objective Function: Volume fraction (topology optimization). */ su2double Total_OFDiscreteness; /*!< \brief Total Objective Function: Discreteness (topology optimization). */ su2double Total_OFCompliance; /*!< \brief Total Objective Function: Compliance (topology optimization). */ + su2double Total_OFStressPenalty; /*!< \brief Total Objective Function: Stress penalty. */ su2double ObjFunc; su2double Global_OFRefGeom; /*!< \brief Global Objective Function (added over time steps): Reference Geometry. */ @@ -546,7 +547,6 @@ class CFEASolver : public CSolver { /*! * \brief Provide the maximum Von Mises Stress for structural analysis. - * \return Value of the maximum Von Mises Stress. */ inline su2double GetTotal_CFEA(void) const final { return Total_CFEA; } @@ -572,10 +572,14 @@ class CFEASolver : public CSolver { /*! * \brief Retrieve the value of the structural compliance objective function - * \return Value of the objective function. */ inline su2double GetTotal_OFCompliance(void) const final { return Total_OFCompliance; } + /*! + * \brief Retrieve the value of the stress penalty objective function + */ + inline su2double GetTotal_OFStressPenalty(void) const final { return Total_OFStressPenalty; } + /*! * \brief Compute the objective function. * \param[in] config - Definition of the problem. @@ -598,6 +602,9 @@ class CFEASolver : public CSolver { case TOPOL_DISCRETENESS: ObjFunc = GetTotal_OFDiscreteness(); break; + case STRESS_PENALTY: + ObjFunc = GetTotal_OFStressPenalty(); + break; } } diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index 2390a546ebfb..a48f946819f5 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -2465,16 +2465,22 @@ class CSolver { inline virtual su2double GetTotal_OFVolFrac() const { return 0; } /*! - * \brief Retrieve the value of the discreteness objective function + * \brief Retrieve the value of the discreteness objective function. */ inline virtual su2double GetTotal_OFDiscreteness() const { return 0; } /*! * \brief A virtual member. - * \return Value of the objective function for the structural compliance. + * \return Value of the compliance objective function. */ inline virtual su2double GetTotal_OFCompliance() const { return 0; } + /*! + * \brief A virtual member. + * \return Value of the stress penalty objective function. + */ + inline virtual su2double GetTotal_OFStressPenalty() const { return 0; } + /*! * \brief A virtual member. * \return Bool that defines whether the solution has an element-based file or not diff --git a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp index 737ea55bbdc1..31c5394577da 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEAElasticity.cpp @@ -243,6 +243,11 @@ void CFEAElasticity::Compute_Dead_Load(CElement *element, const CConfig *config) SetElement_Properties(element, config); /*-----------------------------------------------------------*/ + /*--- Register pre-accumulation inputs, density and reference coords. ---*/ + AD::StartPreacc(); + AD::SetPreaccIn(Rho_s_DL); + element->SetPreaccIn_Coords(false); + unsigned short iGauss, nGauss; unsigned short iNode, iDim, nNode; @@ -286,6 +291,10 @@ void CFEAElasticity::Compute_Dead_Load(CElement *element, const CConfig *config) } + /*--- Register the dead load as preaccumulation output. ---*/ + element->SetPreaccOut_FDL_a(); + AD::EndPreacc(); + } diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp index 44da9a63da88..1f0e6201186e 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp @@ -56,8 +56,6 @@ void CFEALinearElasticity::Compute_Tangent_Matrix(CElement *element, const CConf AD::StartPreacc(); AD::SetPreaccIn(E); AD::SetPreaccIn(Nu); - AD::SetPreaccIn(Rho_s); - AD::SetPreaccIn(Rho_s_DL); element->SetPreaccIn_Coords(); /*--- Recompute Lame parameters as they depend on the material properties ---*/ Compute_Lame_Parameters(); @@ -247,6 +245,18 @@ void CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, const /*--- Set element properties and recompute the constitutive matrix, this is needed for multiple material cases and for correct differentiation ---*/ SetElement_Properties(element, config); + + /*--- Register pre-accumulation inputs ---*/ + /*--- WARNING: Outputs must be registered outside of this method, this allows more + * flexibility in selecting what is captured by AD, capturing the entire stress + * tensor would use more memory than that used by the stress residuals. ---*/ + AD::StartPreacc(); + AD::SetPreaccIn(E); + AD::SetPreaccIn(Nu); + element->SetPreaccIn_Coords(); + /*--- Recompute Lame parameters as they depend on the material properties ---*/ + Compute_Lame_Parameters(); + Compute_Constitutive_Matrix(element, config); /*--- Initialize auxiliary matrices ---*/ diff --git a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp index 135030eacb8c..a02a85294e02 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp @@ -240,8 +240,6 @@ void CFEANonlinearElasticity::Compute_Tangent_Matrix(CElement *element, const CC AD::StartPreacc(); AD::SetPreaccIn(E); AD::SetPreaccIn(Nu); - AD::SetPreaccIn(Rho_s); - AD::SetPreaccIn(Rho_s_DL); if (maxwell_stress) { AD::SetPreaccIn(EFieldMod_Ref); AD::SetPreaccIn(ke_DE); @@ -483,8 +481,6 @@ void CFEANonlinearElasticity::Compute_NodalStress_Term(CElement *element, const AD::StartPreacc(); AD::SetPreaccIn(E); AD::SetPreaccIn(Nu); - AD::SetPreaccIn(Rho_s); - AD::SetPreaccIn(Rho_s_DL); if (maxwell_stress) { AD::SetPreaccIn(EFieldMod_Ref); AD::SetPreaccIn(ke_DE); @@ -756,6 +752,21 @@ void CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *element, co if (maxwell_stress) SetElectric_Properties(element, config); /*-----------------------------------------------------------*/ + /*--- Register pre-accumulation inputs ---*/ + /*--- WARNING: Outputs must be registered outside of this method, this allows more + * flexibility in selecting what is captured by AD, capturing the entire stress + * tensor would use more memory than that used by the stress residuals. ---*/ + AD::StartPreacc(); + AD::SetPreaccIn(E); + AD::SetPreaccIn(Nu); + if (maxwell_stress) { + AD::SetPreaccIn(EFieldMod_Ref); + AD::SetPreaccIn(ke_DE); + } + element->SetPreaccIn_Coords(); + /*--- Recompute Lame parameters as they depend on the material properties ---*/ + Compute_Lame_Parameters(); + su2double Weight, Jac_x; element->ClearStress(); diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index 0f7408b01166..da1460df2b1c 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -136,6 +136,7 @@ void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS SetHistoryOutputValue("REFERENCE_NODE", fea_solver->GetTotal_OFRefNode()); SetHistoryOutputValue("TOPOL_COMPLIANCE", fea_solver->GetTotal_OFCompliance()); + SetHistoryOutputValue("STRESS_PENALTY", fea_solver->GetTotal_OFStressPenalty()); if (config->GetRefGeom()) { SetHistoryOutputValue("REFERENCE_GEOMETRY", fea_solver->GetTotal_OFRefGeom()); } @@ -171,6 +172,7 @@ void CElasticityOutput::SetHistoryOutputFields(CConfig *config){ AddHistoryOutput("REFERENCE_NODE", "RefNode", ScreenOutputFormat::SCIENTIFIC, "STRUCT_COEFF", "", HistoryFieldType::COEFFICIENT); AddHistoryOutput("TOPOL_COMPLIANCE", "TopComp", ScreenOutputFormat::SCIENTIFIC, "STRUCT_COEFF", "", HistoryFieldType::COEFFICIENT); + AddHistoryOutput("STRESS_PENALTY", "StressPen", ScreenOutputFormat::SCIENTIFIC, "STRUCT_COEFF", "", HistoryFieldType::COEFFICIENT); if (config->GetRefGeom()) { AddHistoryOutput("REFERENCE_GEOMETRY", "RefGeom", ScreenOutputFormat::SCIENTIFIC, "STRUCT_COEFF", "", HistoryFieldType::COEFFICIENT); } @@ -178,6 +180,7 @@ void CElasticityOutput::SetHistoryOutputFields(CConfig *config){ AddHistoryOutput("VOLUME_FRACTION", "VolFrac", ScreenOutputFormat::SCIENTIFIC, "STRUCT_COEFF", "", HistoryFieldType::COEFFICIENT); AddHistoryOutput("TOPOL_DISCRETENESS", "TopDisc", ScreenOutputFormat::SCIENTIFIC, "STRUCT_COEFF", "", HistoryFieldType::COEFFICIENT); } + } void CElasticityOutput::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 f9afa46d1651..8e84ceebe29a 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1321,24 +1321,51 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric } -void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { +namespace { +template +su2double vonMissesStress(unsigned short nDim, const T& stress) { + if (nDim == 2) { + su2double Sxx = stress[0], Syy = stress[1], Sxy = stress[2]; - /*--- Never record this method as atm it is not differentiable. ---*/ - const bool wasActive = AD::BeginPassive(); + su2double S1, S2; S1 = S2 = (Sxx+Syy)/2; + su2double tauMax = sqrt(pow((Sxx-Syy)/2, 2) + pow(Sxy,2)); + S1 += tauMax; + S2 -= tauMax; + + return sqrt(S1*S1+S2*S2-2*S1*S2); + } + else { + su2double Sxx = stress[0], Syy = stress[1], Szz = stress[3]; + su2double Sxy = stress[2], Sxz = stress[4], Syz = stress[5]; + + return sqrt(0.5*(pow(Sxx - Syy, 2) + + pow(Syy - Szz, 2) + + pow(Szz - Sxx, 2) + + 6.0*(Sxy*Sxy+Sxz*Sxz+Syz*Syz))); + } +} +} + +void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { const bool prestretch_fem = config->GetPrestretch(); const bool topology_mode = config->GetTopology_Optimization(); - const su2double simp_exponent = config->GetSIMP_Exponent(); + const auto simp_exponent = config->GetSIMP_Exponent(); const unsigned short nStress = (nDim == 2) ? 3 : 6; + const auto stressScale = 1.0 / config->GetAllowedVMStress(); + su2double StressPenalty = 0.0; su2double MaxVonMises_Stress = 0.0; /*--- Start OpenMP parallel region. ---*/ SU2_OMP_PARALLEL { + /*--- Some parts are not recorded, atm only StressPenalty is differentiated to save memory. ---*/ + bool wasActive = AD::BeginPassive(); + /*--- Clear reactions. ---*/ LinSysReact.SetValZero(); @@ -1349,9 +1376,12 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, nodes->SetStress_FEM(iPoint,iStress, 0.0); } } + AD::EndPassive(wasActive); for(auto color : ElemColoring) { + su2double stressPen = 0.0; + /*--- 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) { @@ -1379,7 +1409,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, for (iDim = 0; iDim < nDim; iDim++) { /*--- Compute current coordinate. ---*/ - su2double val_Coord = Get_ValCoord(geometry, indexNode[iNode], iDim); + su2double val_Coord = geometry->nodes->GetCoord(indexNode[iNode],iDim); su2double val_Sol = nodes->GetSolution(indexNode[iNode],iDim) + val_Coord; /*--- If pre-stretched the reference coordinate is stored in the nodes. ---*/ @@ -1406,6 +1436,12 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, numerics[NUM_TERM]->Compute_Averaged_NodalStress(element, config); + /*--- Average the stresses for the element to compute the stress penalty, + * this is the way to make the AD of that function more economical. ---*/ + su2double avgStressElem[6] = {0.0}; + AD::SetPreaccIn(simp_penalty); + const su2double el_weight = simp_penalty / nNodes; + for (iNode = 0; iNode < nNodes; iNode++) { auto iPoint = indexNode[iNode]; @@ -1417,18 +1453,30 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, LinSysReact(iPoint,iVar) += simp_penalty*Ta[iVar]; /*--- Divide the nodal stress by the number of elements that will contribute to this point. ---*/ - su2double weight = simp_penalty / geometry->nodes->GetnElem(iPoint); + su2double pt_weight = simp_penalty / geometry->nodes->GetnElem(iPoint); - for (iStress = 0; iStress < nStress; iStress++) - nodes->AddStress_FEM(iPoint,iStress, weight*element->Get_NodalStress(iNode,iStress)); + for (iStress = 0; iStress < nStress; iStress++) { + auto sigma = element->Get_NodalStress(iNode,iStress); + nodes->AddStress_FEM(iPoint,iStress, pt_weight*sigma); + avgStressElem[iStress] += el_weight*sigma; + } if (LockStrategy) omp_unset_lock(&UpdateLocks[iPoint]); } + su2double elPenalty = pow(max(0.0, vonMissesStress(nDim,avgStressElem)*stressScale - 1.0), 2); + AD::SetPreaccOut(elPenalty); + AD::EndPreacc(); // started above in Compute_Averaged_NodalStress + + stressPen += elPenalty; + } // end iElem loop + SU2_OMP_ATOMIC + StressPenalty += stressPen; } // end color loop + wasActive = AD::BeginPassive(); /*--- Compute the von Misses stress at each point, and the maximum for the domain. ---*/ su2double maxVonMises = 0.0; @@ -1436,51 +1484,24 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, SU2_OMP(for schedule(static,omp_chunk_size) nowait) for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { - /*--- Get the stresses, added up from all the elements that connect to the node. ---*/ - - auto Stress = nodes->GetStress_FEM(iPoint); - su2double VonMises_Stress; - - if (nDim == 2) { - - su2double Sxx = Stress[0], Syy = Stress[1], Sxy = Stress[2]; - - su2double S1, S2; S1 = S2 = (Sxx+Syy)/2; - su2double tauMax = sqrt(pow((Sxx-Syy)/2, 2) + pow(Sxy,2)); - S1 += tauMax; - S2 -= tauMax; + const auto vms = vonMissesStress(nDim, nodes->GetStress_FEM(iPoint)); - VonMises_Stress = sqrt(S1*S1+S2*S2-2*S1*S2); - } - else { - - su2double Sxx = Stress[0], Syy = Stress[1], Szz = Stress[3]; - su2double Sxy = Stress[2], Sxz = Stress[4], Syz = Stress[5]; - - VonMises_Stress = sqrt(0.5*(pow(Sxx - Syy, 2) + - pow(Syy - Szz, 2) + - pow(Szz - Sxx, 2) + - 6.0*(Sxy*Sxy+Sxz*Sxz+Syz*Syz))); - } + nodes->SetVonMises_Stress(iPoint, vms); - nodes->SetVonMises_Stress(iPoint,VonMises_Stress); - - /*--- Update the maximum value of the Von Mises Stress ---*/ - - maxVonMises = max(maxVonMises, VonMises_Stress); + maxVonMises = max(maxVonMises, vms); } SU2_OMP_CRITICAL MaxVonMises_Stress = max(MaxVonMises_Stress, maxVonMises); - } // end SU2_OMP_PARALLEL + AD::EndPassive(wasActive); - su2double tmp = MaxVonMises_Stress; - SU2_MPI::Allreduce(&tmp, &MaxVonMises_Stress, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); + } // end SU2_OMP_PARALLEL /*--- Set the value of the MaxVonMises_Stress as the CFEA coeffient ---*/ + SU2_MPI::Allreduce(&MaxVonMises_Stress, &Total_CFEA, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); - Total_CFEA = MaxVonMises_Stress; - + /*--- Reduce the stress penalty over all ranks ---*/ + SU2_MPI::Allreduce(&StressPenalty, &Total_OFStressPenalty, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); bool outputReactions = false; @@ -1603,8 +1624,6 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, } - AD::EndPassive(wasActive); - } void CFEASolver::Compute_DeadLoad(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { @@ -1928,6 +1947,9 @@ void CFEASolver::Postprocessing(CGeometry *geometry, CConfig *config, CNumerics case VOLUME_FRACTION: Compute_OFVolFrac(geometry, config); break; case TOPOL_DISCRETENESS: Compute_OFVolFrac(geometry, config); break; case TOPOL_COMPLIANCE: Compute_OFCompliance(geometry, config); break; + case STRESS_PENALTY: + Compute_NodalStress(geometry, numerics, config); + break; } return; } From 162265a00536a13a4a9966c52a486dec31b86371 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 19 Jan 2021 17:54:17 +0000 Subject: [PATCH 2/9] fix the AD problem and avoid stopping the preacc outside of what starts it --- SU2_CFD/include/numerics/CNumerics.hpp | 2 +- .../numerics/elasticity/CFEAElasticity.hpp | 28 +++++++++- .../elasticity/CFEALinearElasticity.hpp | 2 +- .../elasticity/CFEANonlinearElasticity.hpp | 2 +- .../elasticity/CFEALinearElasticity.cpp | 19 +++++-- .../elasticity/CFEANonlinearElasticity.cpp | 17 +++++- SU2_CFD/src/solvers/CFEASolver.cpp | 55 ++++--------------- 7 files changed, 71 insertions(+), 54 deletions(-) diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index 368c83cebc1b..29ca32b66b45 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -1495,7 +1495,7 @@ class CNumerics { * \brief A virtual member to compute the averaged nodal stresses * \param[in] element_container - Element structure for the particular element integrated. */ - inline virtual void Compute_Averaged_NodalStress(CElement *element_container, const CConfig* config) { } + inline virtual su2double Compute_Averaged_NodalStress(CElement *element_container, const CConfig* config) { return 0; } /*! * \brief Computes a basis of orthogonal vectors from a supplied vector diff --git a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp index fd96e38a99d9..b1583466cf1c 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEAElasticity.hpp @@ -174,7 +174,33 @@ class CFEAElasticity : public CNumerics { * \param[in,out] element_container - The finite element. * \param[in] config - Definition of the problem. */ - inline void Compute_Averaged_NodalStress(CElement *element_container, const CConfig *config) override { }; + inline su2double Compute_Averaged_NodalStress(CElement *element_container, const CConfig *config) override { return 0; }; + + /*! + * \brief Compute VonMises stress from components Sxx Syy Sxy Szz Sxz Syz. + */ + template + static su2double VonMisesStress(unsigned short nDim, const T& stress) { + if (nDim == 2) { + su2double Sxx = stress[0], Syy = stress[1], Sxy = stress[2]; + + su2double S1, S2; S1 = S2 = (Sxx+Syy)/2; + su2double tauMax = sqrt(pow((Sxx-Syy)/2, 2) + pow(Sxy,2)); + S1 += tauMax; + S2 -= tauMax; + + return sqrt(S1*S1+S2*S2-2*S1*S2); + } + else { + su2double Sxx = stress[0], Syy = stress[1], Szz = stress[3]; + su2double Sxy = stress[2], Sxz = stress[4], Syz = stress[5]; + + return sqrt(0.5*(pow(Sxx - Syy, 2) + + pow(Syy - Szz, 2) + + pow(Szz - Sxx, 2) + + 6.0*(Sxy*Sxy+Sxz*Sxz+Syz*Syz))); + } + } protected: /*! diff --git a/SU2_CFD/include/numerics/elasticity/CFEALinearElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEALinearElasticity.hpp index 7541cc11dd7e..afd58c34e98d 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEALinearElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEALinearElasticity.hpp @@ -72,7 +72,7 @@ class CFEALinearElasticity : public CFEAElasticity { * \param[in,out] element_container - The finite element. * \param[in] config - Definition of the problem. */ - void Compute_Averaged_NodalStress(CElement *element_container, const CConfig *config) final; + su2double Compute_Averaged_NodalStress(CElement *element_container, const CConfig *config) final; private: /*! diff --git a/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp b/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp index 3ef778734aea..8a1f5206711f 100644 --- a/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp +++ b/SU2_CFD/include/numerics/elasticity/CFEANonlinearElasticity.hpp @@ -127,7 +127,7 @@ class CFEANonlinearElasticity : public CFEAElasticity { * \param[in,out] element_container - The finite element. * \param[in] config - Definition of the problem. */ - void Compute_Averaged_NodalStress(CElement *element_container, const CConfig *config) final; + su2double Compute_Averaged_NodalStress(CElement *element_container, const CConfig *config) final; protected: /*! diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp index 1f0e6201186e..84d0f7215e02 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp @@ -232,7 +232,7 @@ void CFEALinearElasticity::Compute_Constitutive_Matrix(CElement *element_contain } -void CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, const CConfig *config) { +su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, const CConfig *config) { unsigned short iVar, jVar; unsigned short iGauss, nGauss; @@ -240,16 +240,13 @@ void CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, const unsigned short iDim, bDim; /*--- Auxiliary vector ---*/ - su2double Strain[6], Stress[6]; + su2double Strain[DIM_STRAIN_3D], Stress[DIM_STRAIN_3D], avgStress[DIM_STRAIN_3D] = {0.0}; /*--- Set element properties and recompute the constitutive matrix, this is needed for multiple material cases and for correct differentiation ---*/ SetElement_Properties(element, config); /*--- Register pre-accumulation inputs ---*/ - /*--- WARNING: Outputs must be registered outside of this method, this allows more - * flexibility in selecting what is captured by AD, capturing the entire stress - * tensor would use more memory than that used by the stress residuals. ---*/ AD::StartPreacc(); AD::SetPreaccIn(E); AD::SetPreaccIn(Nu); @@ -351,6 +348,18 @@ void CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, const } + for (unsigned short iStress = 0; iStress < bDim; ++iStress) + for (iNode = 0; iNode < nNode; iNode++) + avgStress[iStress] += element->Get_NodalStress(iNode, iStress) / nNode; + + auto elStress = VonMisesStress(nDim, avgStress); + + /*--- We only differentiate w.r.t. an avg VM stress for the element as + * considering all nodal stresses would use too much memory. ---*/ + AD::SetPreaccOut(elStress); + AD::EndPreacc(); + + return elStress; } diff --git a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp index a02a85294e02..5c7a338ba7fc 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp @@ -741,7 +741,7 @@ void CFEANonlinearElasticity::Assign_cijkl_D_Mat(void) { } -void CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *element, const CConfig *config) { +su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *element, const CConfig *config) { unsigned short iVar, jVar, kVar; unsigned short iGauss, nGauss; @@ -879,4 +879,19 @@ void CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *element, co } + su2double avgStress[DIM_STRAIN_3D] = {0.0}; + const auto nStress = (nDim == 2) ? DIM_STRAIN_2D : DIM_STRAIN_3D; + + for (unsigned short iStress = 0; iStress < nStress; ++iStress) + for (iNode = 0; iNode < nNode; iNode++) + avgStress[iStress] += element->Get_NodalStress(iNode, iStress) / nNode; + + auto elStress = VonMisesStress(nDim, avgStress); + + /*--- We only differentiate w.r.t. an avg VM stress for the element as + * considering all nodal stresses would use too much memory. ---*/ + AD::SetPreaccOut(elStress); + AD::EndPreacc(); + + return elStress; } diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 8e84ceebe29a..f874260e87f9 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -27,6 +27,7 @@ #include "../../include/solvers/CFEASolver.hpp" #include "../../include/variables/CFEABoundVariable.hpp" +#include "../../include/numerics/elasticity/CFEAElasticity.hpp" #include "../../../Common/include/toolboxes/printing_toolbox.hpp" #include "../../../Common/include/toolboxes/geometry_toolbox.hpp" #include @@ -1321,31 +1322,6 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric } -namespace { -template -su2double vonMissesStress(unsigned short nDim, const T& stress) { - if (nDim == 2) { - su2double Sxx = stress[0], Syy = stress[1], Sxy = stress[2]; - - su2double S1, S2; S1 = S2 = (Sxx+Syy)/2; - su2double tauMax = sqrt(pow((Sxx-Syy)/2, 2) + pow(Sxy,2)); - S1 += tauMax; - S2 -= tauMax; - - return sqrt(S1*S1+S2*S2-2*S1*S2); - } - else { - su2double Sxx = stress[0], Syy = stress[1], Szz = stress[3]; - su2double Sxy = stress[2], Sxz = stress[4], Syz = stress[5]; - - return sqrt(0.5*(pow(Sxx - Syy, 2) + - pow(Syy - Szz, 2) + - pow(Szz - Sxx, 2) + - 6.0*(Sxy*Sxy+Sxz*Sxz+Syz*Syz))); - } -} -} - void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, const CConfig *config) { const bool prestretch_fem = config->GetPrestretch(); @@ -1355,7 +1331,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, const unsigned short nStress = (nDim == 2) ? 3 : 6; - const auto stressScale = 1.0 / config->GetAllowedVMStress(); + const su2double stressScale = 1.0 / config->GetAllowedVMStress(); su2double StressPenalty = 0.0; su2double MaxVonMises_Stress = 0.0; @@ -1434,13 +1410,11 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, /*--- Compute the averaged nodal stresses. ---*/ int NUM_TERM = thread*MAX_TERMS + element_properties[iElem]->GetMat_Mod(); - numerics[NUM_TERM]->Compute_Averaged_NodalStress(element, config); + auto elStress = numerics[NUM_TERM]->Compute_Averaged_NodalStress(element, config); - /*--- Average the stresses for the element to compute the stress penalty, - * this is the way to make the AD of that function more economical. ---*/ - su2double avgStressElem[6] = {0.0}; - AD::SetPreaccIn(simp_penalty); - const su2double el_weight = simp_penalty / nNodes; + stressPen += pow(max(0.0, elStress*simp_penalty*stressScale - 1.0), 2); + + wasActive = AD::BeginPassive(); for (iNode = 0; iNode < nNodes; iNode++) { @@ -1453,22 +1427,15 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, LinSysReact(iPoint,iVar) += simp_penalty*Ta[iVar]; /*--- Divide the nodal stress by the number of elements that will contribute to this point. ---*/ - su2double pt_weight = simp_penalty / geometry->nodes->GetnElem(iPoint); + su2double weight = simp_penalty / geometry->nodes->GetnElem(iPoint); - for (iStress = 0; iStress < nStress; iStress++) { - auto sigma = element->Get_NodalStress(iNode,iStress); - nodes->AddStress_FEM(iPoint,iStress, pt_weight*sigma); - avgStressElem[iStress] += el_weight*sigma; - } + for (iStress = 0; iStress < nStress; iStress++) + nodes->AddStress_FEM(iPoint,iStress, weight*element->Get_NodalStress(iNode,iStress)); if (LockStrategy) omp_unset_lock(&UpdateLocks[iPoint]); } - su2double elPenalty = pow(max(0.0, vonMissesStress(nDim,avgStressElem)*stressScale - 1.0), 2); - AD::SetPreaccOut(elPenalty); - AD::EndPreacc(); // started above in Compute_Averaged_NodalStress - - stressPen += elPenalty; + AD::EndPassive(wasActive); } // end iElem loop SU2_OMP_ATOMIC @@ -1484,7 +1451,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, SU2_OMP(for schedule(static,omp_chunk_size) nowait) for (auto iPoint = 0ul; iPoint < nPointDomain; iPoint++) { - const auto vms = vonMissesStress(nDim, nodes->GetStress_FEM(iPoint)); + const auto vms = CFEAElasticity::VonMisesStress(nDim, nodes->GetStress_FEM(iPoint)); nodes->SetVonMises_Stress(iPoint, vms); From a5e8a7479247a2aadfc987f08cf1d8dcd4d2c524 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 22 Jan 2021 17:32:33 +0000 Subject: [PATCH 3/9] fix COptionDoubleArray, const correctness, and one bug --- Common/include/CConfig.hpp | 205 ++++++-------- Common/include/option_structure.inl | 39 +-- Common/src/CConfig.cpp | 252 ++++++++---------- Common/src/grid_movement/CSurfaceMovement.cpp | 2 +- SU2_CFD/src/solvers/CEulerSolver.cpp | 22 +- SU2_CFD/src/solvers/CFEASolver.cpp | 8 +- SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp | 8 +- SU2_PY/SU2/io/historyMap.py | 17 ++ TestCases/fea_topology/config.cfg | 4 + 9 files changed, 251 insertions(+), 306 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index babf1eb237ad..0ea1ffc0c9e7 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -78,10 +78,7 @@ class CConfig { su2double Fan_Poly_Eff; /*!< \brief Fan polytropic effeciency. */ su2double MinLogResidual; /*!< \brief Minimum value of the log residual. */ su2double EA_ScaleFactor; /*!< \brief Equivalent Area scaling factor */ - su2double* EA_IntLimit; /*!< \brief Integration limits of the Equivalent Area computation */ su2double AdjointLimit; /*!< \brief Adjoint variable limit */ - su2double* Obj_ChainRuleCoeff; /*!< \brief Array defining objective function for adjoint problem based on - chain rule in terms of gradient w.r.t. density, velocity, pressure */ string* ConvField; /*!< \brief Field used for convergence check.*/ string* WndConvField; /*!< \brief Function where to apply the windowed convergence criteria for the time average of the unsteady (single zone) flow problem. */ @@ -93,10 +90,6 @@ class CConfig { bool Wnd_Cauchy_Crit; /*!< \brief True => Cauchy criterion is used for time average objective function in unsteady flows. */ bool MG_AdjointFlow; /*!< \brief MG with the adjoint flow problem */ - su2double* SubsonicEngine_Cyl; /*!< \brief Coordinates of the box subsonic region */ - su2double* SubsonicEngine_Values; /*!< \brief Values of the box subsonic region */ - su2double* Hold_GridFixed_Coord; /*!< \brief Coordinates of the box to hold fixed the nbumerical grid */ - su2double *DistortionRack; su2double *PressureLimits, *DensityLimits, *TemperatureLimits; /*!< \brief Limits for the primitive variables */ @@ -106,7 +99,6 @@ class CConfig { unsigned short ConvCriteria; /*!< \brief Kind of convergence criteria. */ unsigned short nFFD_Iter; /*!< \brief Iteration for the point inversion problem. */ unsigned short FFD_Blending; /*!< \brief Kind of FFD Blending function. */ - su2double* FFD_BSpline_Order; /*!< \brief BSpline order in i,j,k direction. */ su2double FFD_Tol; /*!< \brief Tolerance in the point inversion problem. */ bool FFD_IntPrev; /*!< \brief Enables self-intersection prevention procedure within the FFD box. */ unsigned short FFD_IntPrev_MaxIter; /*!< \brief Amount of iterations for FFD box self-intersection prevention procedure. */ @@ -470,7 +462,6 @@ class CConfig { *MG_PostSmooth, /*!< \brief Multigrid Post smoothing. */ *MG_CorrecSmooth; /*!< \brief Multigrid Jacobi implicit smoothing of the correction. */ su2double *LocationStations; /*!< \brief Airfoil sections in wing slicing subroutine. */ - su2double *NacelleLocation; /*!< \brief Definition of the nacelle location. */ unsigned short Kind_Solver, /*!< \brief Kind of solver Euler, NS, Continuous adjoint, etc. */ Kind_MZSolver, /*!< \brief Kind of multizone solver. */ @@ -592,13 +583,8 @@ class CConfig { su2double AdjTurb_Linear_Error; /*!< \brief Min error of the turbulent adjoint linear solver for the implicit formulation. */ su2double EntropyFix_Coeff; /*!< \brief Entropy fix coefficient. */ unsigned short AdjTurb_Linear_Iter; /*!< \brief Min error of the turbulent adjoint linear solver for the implicit formulation. */ - su2double *Stations_Bounds; /*!< \brief Airfoil section limit. */ unsigned short nLocationStations, /*!< \brief Number of section cuts to make when outputting mesh and cp . */ nWingStations; /*!< \brief Number of section cuts to make when calculating internal volume. */ - su2double* Kappa_Flow, /*!< \brief Numerical dissipation coefficients for the flow equations. */ - *Kappa_AdjFlow, /*!< \brief Numerical dissipation coefficients for the adjoint flow equations. */ - *Kappa_Heat; /*!< \brief Numerical dissipation coefficients for the (fvm) heat equation. */ - su2double* FFD_Axis; /*!< \brief Numerical dissipation coefficients for the adjoint equations. */ su2double Kappa_1st_AdjFlow, /*!< \brief Lax 1st order dissipation coefficient for adjoint flow equations (coarse multigrid levels). */ Kappa_2nd_AdjFlow, /*!< \brief JST 2nd order dissipation coefficient for adjoint flow equations. */ Kappa_4th_AdjFlow, /*!< \brief JST 4th order dissipation coefficient for adjoint flow equations. */ @@ -747,7 +733,6 @@ class CConfig { *CFL_AdaptParam, /*!< \brief Information about the CFL ramp. */ *RelaxFactor_Giles, /*!< \brief Information about the under relaxation factor for Giles BC. */ *CFL, /*!< \brief CFL number. */ - *HTP_Axis, /*!< \brief Location of the HTP axis. */ DomainVolume; /*!< \brief Volume of the computational grid. */ unsigned short nRefOriginMoment_X, /*!< \brief Number of X-coordinate moment computation origins. */ @@ -755,8 +740,6 @@ class CConfig { nRefOriginMoment_Z; /*!< \brief Number of Z-coordinate moment computation origins. */ unsigned short nMesh_Box_Size; short *Mesh_Box_Size; /*!< \brief Array containing the number of grid points in the x-, y-, and z-directions for the analytic RECTANGLE and BOX grid formats. */ - su2double* Mesh_Box_Length; /*!< \brief Array containing the length in the x-, y-, and z-directions for the analytic RECTANGLE and BOX grid formats. */ - su2double* Mesh_Box_Offset; /*!< \brief Array containing the offset from 0.0 in the x-, y-, and z-directions for the analytic RECTANGLE and BOX grid formats. */ string Mesh_FileName, /*!< \brief Mesh input file. */ Mesh_Out_FileName, /*!< \brief Mesh output file. */ Solution_FileName, /*!< \brief Flow solution input file. */ @@ -799,7 +782,6 @@ class CConfig { Inc_Velocity_Ref, /*!< \brief Reference velocity for custom incompressible non-dim. */ Inc_Temperature_Ref, /*!< \brief Reference temperature for custom incompressible non-dim. */ Inc_Density_Init, /*!< \brief Initial density for incompressible flows. */ - *Inc_Velocity_Init, /*!< \brief Initial velocity vector for incompressible flows. */ Inc_Temperature_Init, /*!< \brief Initial temperature for incompressible flows w/ heat transfer. */ Heat_Flux_Ref, /*!< \brief Reference heat flux for non-dim. */ Gas_Constant_Ref, /*!< \brief Reference specific gas constant. */ @@ -817,9 +799,6 @@ class CConfig { Mu_Temperature_RefND, /*!< \brief Non-dimensional reference temperature for Sutherland model. */ Mu_S, /*!< \brief Reference S for Sutherland model. */ Mu_SND; /*!< \brief Non-dimensional reference S for Sutherland model. */ - su2double* CpPolyCoefficients; /*!< \brief Definition of the temperature polynomial coefficients for specific heat Cp. */ - su2double* MuPolyCoefficients; /*!< \brief Definition of the temperature polynomial coefficients for viscosity. */ - su2double* KtPolyCoefficients; /*!< \brief Definition of the temperature polynomial coefficients for thermal conductivity. */ array CpPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for specific heat Cp. */ array MuPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for viscosity. */ array KtPolyCoefficientsND{{0.0}}; /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for thermal conductivity. */ @@ -827,7 +806,6 @@ class CConfig { Thermal_Diffusivity_Solid, /*!< \brief Thermal diffusivity in solids. */ Temperature_Freestream_Solid, /*!< \brief Temperature in solids at freestream conditions. */ Density_Solid, /*!< \brief Total density in solids. */ - *Velocity_FreeStream, /*!< \brief Free-stream velocity vector of the fluid. */ Energy_FreeStream, /*!< \brief Free-stream total energy of the fluid. */ ModVel_FreeStream, /*!< \brief Magnitude of the free-stream velocity of the fluid. */ ModVel_FreeStreamND, /*!< \brief Non-dimensional magnitude of the free-stream velocity of the fluid. */ @@ -888,20 +866,19 @@ class CConfig { su2double AitkenDynMinInit; /*!< \brief Aitken's minimum dynamic relaxation factor for the first iteration */ bool RampAndRelease; /*!< \brief option for ramp load and release */ bool Sine_Load; /*!< \brief option for sine load */ - su2double *SineLoad_Coeff; /*!< \brief Stores the load coefficient */ su2double Thermal_Diffusivity; /*!< \brief Thermal diffusivity used in the heat solver. */ su2double Cyclic_Pitch, /*!< \brief Cyclic pitch for rotorcraft simulations. */ Collective_Pitch; /*!< \brief Collective pitch for rotorcraft simulations. */ su2double Mach_Motion; /*!< \brief Mach number based on mesh velocity and freestream quantities. */ - su2double *Motion_Origin, /*!< \brief Mesh motion origin. */ - *Translation_Rate, /*!< \brief Translational velocity of the mesh. */ - *Rotation_Rate, /*!< \brief Angular velocity of the mesh . */ - *Pitching_Omega, /*!< \brief Angular frequency of the mesh pitching. */ - *Pitching_Ampl, /*!< \brief Pitching amplitude. */ - *Pitching_Phase, /*!< \brief Pitching phase offset. */ - *Plunging_Omega, /*!< \brief Angular frequency of the mesh plunging. */ - *Plunging_Ampl; /*!< \brief Plunging amplitude. */ + su2double Motion_Origin[3] = {0.0}, /*!< \brief Mesh motion origin. */ + Translation_Rate[3] = {0.0}, /*!< \brief Translational velocity of the mesh. */ + Rotation_Rate[3] = {0.0}, /*!< \brief Angular velocity of the mesh . */ + Pitching_Omega[3] = {0.0}, /*!< \brief Angular frequency of the mesh pitching. */ + Pitching_Ampl[3] = {0.0}, /*!< \brief Pitching amplitude. */ + Pitching_Phase[3] = {0.0}, /*!< \brief Pitching phase offset. */ + Plunging_Omega[3] = {0.0}, /*!< \brief Angular frequency of the mesh plunging. */ + Plunging_Ampl[3] = {0.0}; /*!< \brief Plunging amplitude. */ su2double *MarkerMotion_Origin, /*!< \brief Mesh motion origin of marker. */ *MarkerTranslation_Rate, /*!< \brief Translational velocity of marker. */ *MarkerRotation_Rate, /*!< \brief Angular velocity of marker. */ @@ -972,7 +949,6 @@ class CConfig { unsigned short Dynamic_LoadTransfer; /*!< \brief Method for dynamic load transferring. */ bool IncrementalLoad; /*!< \brief Apply the load in increments (for nonlinear structural analysis). */ unsigned long IncLoad_Nincrements; /*!< \brief Number of increments. */ - su2double *IncLoad_Criteria; /*!< \brief Criteria for the application of incremental loading. */ su2double Ramp_Time; /*!< \brief Time until the maximum load is applied. */ bool Predictor, /*!< \brief Determines whether a predictor step is used. */ Relaxation; /*!< \brief Determines whether a relaxation step is used. */ @@ -990,15 +966,15 @@ class CConfig { bool FEAAdvancedMode; /*!< \brief Determine if advanced features are used from the element-based FEA analysis (experimental). */ 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. */ - AllowedVMStress; /*!< \brief Maximum stress for the stress penalty objective function. */ + DV_Penalty; /*!< \brief Penalty weight to add a constraint to the total amount of stiffness. */ + array StressPenaltyParam = {{1.0, 20.0}}; /*!< \brief Allowed stress and KS aggregation exponent. */ 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. */ su2double ParMETIS_tolerance; /*!< \brief Load balancing tolerance for ParMETIS. */ long ParMETIS_pointWgt; /*!< \brief Load balancing weight given to points. */ long ParMETIS_edgeWgt; /*!< \brief Load balancing weight given to edges. */ unsigned short DirectDiff; /*!< \brief Direct Differentation mode. */ - bool DiscreteAdjoint; /*!< \brief AD-based discrete adjoint mode. */ + bool DiscreteAdjoint; /*!< \brief AD-based discrete adjoint mode. */ su2double Const_DES; /*!< \brief Detached Eddy Simulation Constant. */ unsigned short Kind_WindowFct; /*!< \brief Type of window (weight) function for objective functional. */ unsigned short Kind_HybridRANSLES; /*!< \brief Kind of Hybrid RANS/LES. */ @@ -1012,19 +988,14 @@ class CConfig { bool SpatialFourier; /*!< \brief option for computing the fourier transforms for subsonic non-reflecting BC. */ bool RampRotatingFrame; /*!< \brief option for ramping up or down the Rotating Frame values */ bool RampOutletPressure; /*!< \brief option for ramping up or down the outlet pressure */ - su2double *Mixedout_Coeff; /*!< \brief coefficient for the */ - su2double *RampRotatingFrame_Coeff; /*!< \brief coefficient for Rotating frame ramp */ - su2double *RampOutletPressure_Coeff; /*!< \brief coefficient for outlet pressure ramp */ su2double AverageMachLimit; /*!< \brief option for turbulent mixingplane */ su2double FinalRotation_Rate_Z; /*!< \brief Final rotation rate Z if Ramp rotating frame is activated. */ su2double FinalOutletPressure; /*!< \brief Final outlet pressure if Ramp outlet pressure is activated. */ su2double MonitorOutletPressure; /*!< \brief Monitor outlet pressure if Ramp outlet pressure is activated. */ - array default_cp_polycoeffs{{0.0}}; /*!< \brief Array for specific heat polynomial coefficients. */ - array default_mu_polycoeffs{{0.0}}; /*!< \brief Array for viscosity polynomial coefficients. */ - array default_kt_polycoeffs{{0.0}}; /*!< \brief Array for thermal conductivity polynomial coefficients. */ - su2double *ExtraRelFacGiles; /*!< \brief coefficient for extra relaxation factor for Giles BC*/ + array cp_polycoeffs{{0.0}}; /*!< \brief Array for specific heat polynomial coefficients. */ + array mu_polycoeffs{{0.0}}; /*!< \brief Array for viscosity polynomial coefficients. */ + array kt_polycoeffs{{0.0}}; /*!< \brief Array for thermal conductivity polynomial coefficients. */ bool Body_Force; /*!< \brief Flag to know if a body force is included in the formulation. */ - su2double *Body_Force_Vector; /*!< \brief Values of the prescribed body force vector. */ su2double *FreeStreamTurboNormal; /*!< \brief Direction to initialize the flow in turbomachinery computation */ su2double Restart_Bandwidth_Agg; /*!< \brief The aggregate of the bandwidth for writing binary restarts (to be averaged later). */ su2double Max_Vel2; /*!< \brief The maximum velocity^2 in the domain for the incompressible preconditioner. */ @@ -1041,11 +1012,9 @@ class CConfig { *top_optim_filter_radius; /*!< \brief Radius of the filter(s) used on the design density for topology optimization. */ unsigned short top_optim_proj_type; /*!< \brief The projection function used in topology optimization. */ su2double top_optim_proj_param; /*!< \brief The value of the parameter for the projection function. */ - bool HeatSource; /*!< \brief Flag to know if there is a volumetric heat source on the flow. */ - su2double ValHeatSource; /*!< \brief Value of the volumetric heat source on the flow (W/m3). */ - su2double Heat_Source_Rot_Z; /*!< \brief Rotation of the volumetric heat source on the Z axis. */ - su2double *Heat_Source_Center, /*!< \brief Position of the center of the heat source. */ - *Heat_Source_Axes; /*!< \brief Principal axes (x, y, z) of the ellipsoid containing the heat source. */ + bool HeatSource; /*!< \brief Flag to know if there is a volumetric heat source on the flow. */ + su2double ValHeatSource; /*!< \brief Value of the volumetric heat source on the flow (W/m3). */ + su2double Heat_Source_Rot_Z; /*!< \brief Rotation of the volumetric heat source on the Z axis. */ unsigned short Kind_Radiation; /*!< \brief Kind of radiation model used. */ unsigned short Kind_P1_Init; /*!< \brief Kind of initialization used in the P1 model. */ su2double Absorption_Coeff, /*!< \brief Absorption coefficient of the medium (radiation). */ @@ -1057,33 +1026,33 @@ class CConfig { su2double CFL_Rad; /*!< \brief CFL Number for the radiation solver. */ array default_cfl_adapt; /*!< \brief Default CFL adapt param array for the COption class. */ - su2double default_vel_inf[3], /*!< \brief Default freestream velocity array for the COption class. */ - default_eng_cyl[7], /*!< \brief Default engine box array for the COption class. */ - default_eng_val[5], /*!< \brief Default engine box array values for the COption class. */ - default_jst_coeff[2], /*!< \brief Default artificial dissipation (flow) array for the COption class. */ - default_ffd_coeff[3], /*!< \brief Default artificial dissipation (flow) array for the COption class. */ - default_mixedout_coeff[3], /*!< \brief Default default mixedout algorithm coefficients for the COption class. */ - default_rampRotFrame_coeff[3], /*!< \brief Default ramp rotating frame coefficients for the COption class. */ - default_rampOutPres_coeff[3], /*!< \brief Default ramp outlet pressure coefficients for the COption class. */ - default_jst_adj_coeff[2], /*!< \brief Default artificial dissipation (adjoint) array for the COption class. */ - default_ad_coeff_heat[2], /*!< \brief Default artificial dissipation (heat) array for the COption class. */ - default_obj_coeff[5], /*!< \brief Default objective array for the COption class. */ - default_mesh_box_length[3], /*!< \brief Default mesh box length for the COption class. */ - default_mesh_box_offset[3], /*!< \brief Default mesh box offset for the COption class. */ - default_geo_loc[2], /*!< \brief Default SU2_GEO section locations array for the COption class. */ - default_distortion[2], /*!< \brief Default SU2_GEO section locations array for the COption class. */ - default_ea_lim[3], /*!< \brief Default equivalent area limit array for the COption class. */ - default_grid_fix[6], /*!< \brief Default fixed grid (non-deforming region) array for the COption class. */ - default_htp_axis[2], /*!< \brief Default HTP axis for the COption class. */ - default_ffd_axis[3], /*!< \brief Default FFD axis for the COption class. */ - default_inc_crit[3], /*!< \brief Default incremental criteria array for the COption class. */ - default_extrarelfac[2], /*!< \brief Default extra relaxation factor for Giles BC in the COption class. */ - default_sineload_coeff[3], /*!< \brief Default values for a sine load. */ - default_body_force[3], /*!< \brief Default body force vector for the COption class. */ - default_nacelle_location[5], /*!< \brief Location of the nacelle. */ - default_hs_axes[3], /*!< \brief Default principal axes (x, y, z) of the ellipsoid containing the heat source. */ - default_hs_center[3], /*!< \brief Default position of the center of the heat source. */ - default_roughness[1]; + su2double vel_init[3], /*!< \brief initial velocity array for the COption class. */ + vel_inf[3], /*!< \brief freestream velocity array for the COption class. */ + eng_cyl[7], /*!< \brief engine box array for the COption class. */ + eng_val[5], /*!< \brief engine box array values for the COption class. */ + jst_coeff[2], /*!< \brief artificial dissipation (flow) array for the COption class. */ + ffd_coeff[3], /*!< \brief artificial dissipation (flow) array for the COption class. */ + mixedout_coeff[3], /*!< \brief default mixedout algorithm coefficients for the COption class. */ + rampRotFrame_coeff[3], /*!< \brief ramp rotating frame coefficients for the COption class. */ + rampOutPres_coeff[3], /*!< \brief ramp outlet pressure coefficients for the COption class. */ + jst_adj_coeff[2], /*!< \brief artificial dissipation (adjoint) array for the COption class. */ + ad_coeff_heat[2], /*!< \brief artificial dissipation (heat) array for the COption class. */ + obj_coeff[5], /*!< \brief objective array for the COption class. */ + mesh_box_length[3], /*!< \brief mesh box length for the COption class. */ + mesh_box_offset[3], /*!< \brief mesh box offset for the COption class. */ + geo_loc[2], /*!< \brief SU2_GEO section locations array for the COption class. */ + distortion[2], /*!< \brief SU2_GEO section locations array for the COption class. */ + ea_lim[3], /*!< \brief equivalent area limit array for the COption class. */ + grid_fix[6], /*!< \brief fixed grid (non-deforming region) array for the COption class. */ + htp_axis[2], /*!< \brief HTP axis for the COption class. */ + ffd_axis[3], /*!< \brief FFD axis for the COption class. */ + inc_crit[3], /*!< \brief incremental criteria array for the COption class. */ + extrarelfac[2], /*!< \brief extra relaxation factor for Giles BC in the COption class. */ + sineload_coeff[3], /*!< \brief values for a sine load. */ + body_force[3], /*!< \brief body force vector for the COption class. */ + nacelle_location[5], /*!< \brief Location of the nacelle. */ + hs_axes[3], /*!< \brief principal axes (x, y, z) of the ellipsoid containing the heat source. */ + hs_center[3]; /*!< \brief position of the center of the heat source. */ unsigned short Riemann_Solver_FEM; /*!< \brief Riemann solver chosen for the DG method. */ su2double Quadrature_Factor_Straight; /*!< \brief Factor applied during quadrature of elements with a constant Jacobian. */ @@ -1219,7 +1188,7 @@ class CConfig { template void addEnumListOption(const string name, unsigned short & input_size, unsigned short * & option_field, const map & enum_map); - void addDoubleArrayOption(const string name, const int size, su2double * & option_field, su2double * default_value); + void addDoubleArrayOption(const string name, const int size, su2double* option_field); void addDoubleListOption(const string name, unsigned short & size, su2double * & option_field); @@ -1450,7 +1419,7 @@ class CConfig { * \param[in] index - 0 means x_min, and 1 means x_max. * \return Integration limits for the equivalent area computation. */ - su2double GetEA_IntLimit(unsigned short index) const { return EA_IntLimit[index]; } + su2double GetEA_IntLimit(unsigned short index) const { return ea_lim[index]; } /*! * \brief Get the integration limits for the equivalent area computation. @@ -1469,25 +1438,25 @@ class CConfig { * \brief Get the coordinates where of the box where the grid is going to be deformed. * \return Coordinates where of the box where the grid is going to be deformed. */ - const su2double *GetHold_GridFixed_Coord(void) const { return Hold_GridFixed_Coord; } + const su2double *GetHold_GridFixed_Coord(void) const { return grid_fix; } /*! * \brief Get the values of subsonic engine. * \return Values of subsonic engine. */ - su2double *GetSubsonicEngine_Values(void) { return SubsonicEngine_Values; } + const su2double *GetSubsonicEngine_Values(void) const { return eng_val; } /*! * \brief Get the cycle of a subsonic engine. * \return Cyl of a subsonic engine. */ - su2double *GetSubsonicEngine_Cyl(void) { return SubsonicEngine_Cyl; } + const su2double *GetSubsonicEngine_Cyl(void) const { return eng_cyl; } /*! * \brief Get the distortion rack. * \return Distortion rack. */ - su2double *GetDistortionRack(void) { return DistortionRack; } + const su2double *GetDistortionRack(void) const { return distortion; } /*! * \brief Get the power of the dual volume in the grid adaptation sensor. @@ -1548,19 +1517,19 @@ class CConfig { * \brief Get the values of the CFL adapation. * \return Value of CFL adapation */ - su2double GetHTP_Axis(unsigned short val_index) const { return HTP_Axis[val_index]; } + su2double GetHTP_Axis(unsigned short val_index) const { return htp_axis[val_index]; } /*! * \brief Get the value of the limits for the sections. * \return Value of the limits for the sections. */ - su2double GetStations_Bounds(unsigned short val_var) const { return Stations_Bounds[val_var]; } + su2double GetStations_Bounds(unsigned short val_var) const { return geo_loc[val_var]; } /*! * \brief Get the value of the vector that connects the cartesian axis with a sherical or cylindrical one. * \return Coordinate of the Axis. */ - su2double GetFFD_Axis(unsigned short val_var) const { return FFD_Axis[val_var]; } + su2double GetFFD_Axis(unsigned short val_var) const { return ffd_axis[val_var]; } /*! * \brief Get the value of the bulk modulus. @@ -1830,8 +1799,8 @@ class CConfig { * \brief Get the vector of the dimensionalized freestream velocity. * \return Dimensionalized freestream velocity vector. */ - su2double* GetVelocity_FreeStream(void) { return Velocity_FreeStream; } - const su2double* GetVelocity_FreeStream(void) const { return Velocity_FreeStream; } + su2double* GetVelocity_FreeStream(void) { return vel_inf; } + const su2double* GetVelocity_FreeStream(void) const { return vel_inf; } /*! * \brief Get the value of the non-dimensionalized freestream temperature. @@ -2020,7 +1989,7 @@ class CConfig { * \brief Get the value of the initial velocity for incompressible flows. * \return Initial velocity for incompressible flows. */ - su2double* GetInc_Velocity_Init(void) { return Inc_Velocity_Init; } + const su2double* GetInc_Velocity_Init(void) const { return vel_init; } /*! * \brief Get the value of the initial temperature for incompressible flows. @@ -2501,7 +2470,7 @@ class CConfig { * \param[in] val_velocity_freestream - Value of the free-stream velocity component. * \param[in] val_dim - Value of the current dimension. */ - void SetVelocity_FreeStream(su2double val_velocity_freestream, unsigned short val_dim) { Velocity_FreeStream[val_dim] = val_velocity_freestream; } + void SetVelocity_FreeStream(su2double val_velocity_freestream, unsigned short val_dim) { vel_inf[val_dim] = val_velocity_freestream; } /*! * \brief Set the Froude number for free surface problems. @@ -2776,7 +2745,7 @@ class CConfig { * \brief Get the kind BSpline Order in i,j,k direction. * \return The kind BSpline Order in i,j,k direction. */ - su2double* GetFFD_BSplineOrder() { return FFD_BSpline_Order;} + const su2double* GetFFD_BSplineOrder() const { return ffd_coeff;} /*! * \brief Get the number of Runge-Kutta steps. @@ -2806,7 +2775,7 @@ class CConfig { * \brief Get the location of the time DOFs for ADER-DG on the interval [-1..1]. * \return The location of the time DOFs used in ADER-DG. */ - su2double *GetTimeDOFsADER_DG(void) { return TimeDOFsADER_DG; } + const su2double *GetTimeDOFsADER_DG(void) const { return TimeDOFsADER_DG; } /*! * \brief Get the number time integration points for ADER-DG. @@ -2818,13 +2787,13 @@ class CConfig { * \brief Get the location of the time integration points for ADER-DG on the interval [-1..1]. * \return The location of the time integration points used in ADER-DG. */ - su2double *GetTimeIntegrationADER_DG(void) { return TimeIntegrationADER_DG; } + const su2double *GetTimeIntegrationADER_DG(void) const { return TimeIntegrationADER_DG; } /*! * \brief Get the weights of the time integration points for ADER-DG. * \return The weights of the time integration points used in ADER-DG. */ - su2double *GetWeightsIntegrationADER_DG(void) { return WeightsIntegrationADER_DG; } + const su2double *GetWeightsIntegrationADER_DG(void) const { return WeightsIntegrationADER_DG; } /*! * \brief Get the total number of boundary markers including send/receive domains. @@ -3469,7 +3438,7 @@ class CConfig { * \param[in] val_index - Index of the section. * \return Coordinate of the nacelle location. */ - su2double GetNacelleLocation(unsigned short val_index) const { return NacelleLocation[val_index]; } + su2double GetNacelleLocation(unsigned short val_index) const { return nacelle_location[val_index]; } /*! * \brief Get the number of pre-smoothings in a multigrid strategy. @@ -3803,7 +3772,7 @@ class CConfig { * \param[in] val_index - Index of the array with all polynomial coefficients. * \return Temperature polynomial coefficient for specific heat Cp. */ - su2double GetCp_PolyCoeff(unsigned short val_index) const { return CpPolyCoefficients[val_index]; } + su2double GetCp_PolyCoeff(unsigned short val_index) const { return cp_polycoeffs[val_index]; } /*! * \brief Get the temperature polynomial coefficient for specific heat Cp. @@ -3817,7 +3786,7 @@ class CConfig { * \param[in] val_index - Index of the array with all polynomial coefficients. * \return Temperature polynomial coefficient for viscosity. */ - su2double GetMu_PolyCoeff(unsigned short val_index) const { return MuPolyCoefficients[val_index]; } + su2double GetMu_PolyCoeff(unsigned short val_index) const { return mu_polycoeffs[val_index]; } /*! * \brief Get the temperature polynomial coefficient for viscosity. @@ -3837,7 +3806,7 @@ class CConfig { * \param[in] val_index - Index of the array with all polynomial coefficients. * \return Temperature polynomial coefficient for thermal conductivity. */ - su2double GetKt_PolyCoeff(unsigned short val_index) const { return KtPolyCoefficients[val_index]; } + su2double GetKt_PolyCoeff(unsigned short val_index) const { return kt_polycoeffs[val_index]; } /*! * \brief Get the temperature polynomial coefficient for thermal conductivity. @@ -4769,7 +4738,7 @@ class CConfig { * \brief Get coeff for Rotating Frame Ramp. * \return coeff Ramp Rotating Frame. */ - su2double GetRampRotatingFrame_Coeff(unsigned short iCoeff) const { return RampRotatingFrame_Coeff[iCoeff];} + su2double GetRampRotatingFrame_Coeff(unsigned short iCoeff) const { return rampRotFrame_coeff[iCoeff];} /*! * \brief Get Rotating Frame Ramp option. @@ -4781,7 +4750,7 @@ class CConfig { * \brief Get coeff for Outlet Pressure Ramp. * \return coeff Ramp Outlet Pressure. */ - su2double GetRampOutletPressure_Coeff(unsigned short iCoeff) const { return RampOutletPressure_Coeff[iCoeff];} + su2double GetRampOutletPressure_Coeff(unsigned short iCoeff) const { return rampOutPres_coeff[iCoeff];} /*! * \brief Get final Outlet Pressure value for the ramp. @@ -4810,13 +4779,13 @@ class CConfig { * \brief Get mixedout coefficients. * \return mixedout coefficient. */ - su2double GetMixedout_Coeff(unsigned short iCoeff) const { return Mixedout_Coeff[iCoeff];} + su2double GetMixedout_Coeff(unsigned short iCoeff) const { return mixedout_coeff[iCoeff];} /*! * \brief Get extra relaxation factor coefficients for the Giels BC. * \return mixedout coefficient. */ - su2double GetExtraRelFacGiles(unsigned short iCoeff) const { return ExtraRelFacGiles[iCoeff];} + su2double GetExtraRelFacGiles(unsigned short iCoeff) const { return extrarelfac[iCoeff];} /*! * \brief Get mach limit for average massflow-based procedure . @@ -5059,7 +5028,7 @@ class CConfig { * calculated using the area averaged outlet values of density, velocity, and pressure. * Gradients are w.r.t density, velocity[3], and pressure. when 2D gradient w.r.t. 3rd component of velocity set to 0. */ - su2double GetCoeff_ObjChainRule(unsigned short iVar) const { return Obj_ChainRuleCoeff[iVar]; } + su2double GetCoeff_ObjChainRule(unsigned short iVar) const { return obj_coeff[iVar]; } /*! * \brief Get the kind of sensitivity smoothing technique. @@ -5692,7 +5661,7 @@ class CConfig { * \brief Get the Harmonic Balance frequency pointer. * \return Harmonic Balance Frequency pointer. */ - su2double* GetOmega_HB(void) { return Omega_HB; } + const su2double* GetOmega_HB(void) const { return Omega_HB; } /*! * \brief Get if harmonic balance source term is to be preconditioned @@ -5765,7 +5734,7 @@ class CConfig { * \brief Get a pointer to the body force vector. * \return A pointer to the body force vector. */ - const su2double* GetBody_Force_Vector(void) const { return Body_Force_Vector; } + const su2double* GetBody_Force_Vector(void) const { return body_force; } /*! * \brief Get information about the volumetric heat source. @@ -5795,7 +5764,7 @@ class CConfig { * \brief Get the position of the center of the volumetric heat source. * \return Pointer to the center of the ellipsoid that introduces a volumetric heat source. */ - inline const su2double* GetHeatSource_Center(void) const {return Heat_Source_Center;} + inline const su2double* GetHeatSource_Center(void) const {return hs_center;} /*! * \brief Set the position of the center of the volumetric heat source. @@ -5804,14 +5773,14 @@ class CConfig { * \param[in] z_cent = Z position of the center of the volumetric heat source. */ inline void SetHeatSource_Center(su2double x_cent, su2double y_cent, su2double z_cent) { - Heat_Source_Center[0] = x_cent; Heat_Source_Center[1] = y_cent; Heat_Source_Center[2] = z_cent; + hs_center[0] = x_cent; hs_center[1] = y_cent; hs_center[2] = z_cent; } /*! * \brief Get the radius of the ellipsoid that introduces a volumetric heat source. * \return Pointer to the radii (x, y, z) of the ellipsoid that introduces a volumetric heat source. */ - inline const su2double* GetHeatSource_Axes(void) const {return Heat_Source_Axes;} + inline const su2double* GetHeatSource_Axes(void) const {return hs_axes;} /*! * \brief Get information about the rotational frame. @@ -6397,7 +6366,7 @@ class CConfig { * \param[in] val_index - Index corresponding to the inlet boundary. * \return The inlet velocity vector. */ - su2double* GetInlet_Velocity(string val_index); + const su2double* GetInlet_Velocity(string val_index) const; /*! * \brief Get the mass fraction vector at a supersonic inlet boundary. @@ -6473,7 +6442,7 @@ class CConfig { * \param[in] val_marker - Index corresponding to the Riemann boundary. * \return The Flowdir */ - su2double* GetRiemann_FlowDir(string val_marker); + const su2double* GetRiemann_FlowDir(string val_marker) const; /*! * \brief Get Kind Data of Riemann boundary. @@ -6501,7 +6470,7 @@ class CConfig { * \param[in] val_marker - Index corresponding to the Giles BC. * \return The Flowdir */ - su2double* GetGiles_FlowDir(string val_marker); + const su2double* GetGiles_FlowDir(string val_marker) const; /*! * \brief Get Kind Data for the Giles BC. @@ -6629,7 +6598,7 @@ class CConfig { * \param[in] val_marker - String of the viscous wall marker. * \return Pointer to the integer info for the given marker. */ - unsigned short* GetWallFunction_IntInfo(string val_marker); + const unsigned short* GetWallFunction_IntInfo(string val_marker) const; /*! * \brief Get the additional double info for the wall function treatment @@ -6637,7 +6606,7 @@ class CConfig { * \param[in] val_marker - String of the viscous wall marker. * \return Pointer to the double info for the given marker. */ - su2double* GetWallFunction_DoubleInfo(string val_marker); + const su2double* GetWallFunction_DoubleInfo(string val_marker) const; /*! * \brief Get the type of wall and roughness height on a wall boundary (Heatflux or Isothermal). @@ -8291,7 +8260,7 @@ class CConfig { * * \brief Set freestream turbonormal for initializing solution. */ - su2double* GetFreeStreamTurboNormal(void) { return FreeStreamTurboNormal; } + const su2double* GetFreeStreamTurboNormal(void) const { return FreeStreamTurboNormal; } /*! * @@ -8511,7 +8480,7 @@ class CConfig { * \param[in] val_index - Index corresponding to the load boundary. * \return The pointer to the sine load values. */ - const su2double* GetLoad_Sine(void) const { 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 @@ -8534,9 +8503,9 @@ class CConfig { su2double GetTotalDV_Penalty(void) const { return DV_Penalty; } /*! - * \brief Get the maximum allowed VM stress for the stress penalty objective function. + * \brief Get the maximum allowed VM stress and KS exponent for the stress penalty objective function. */ - su2double GetAllowedVMStress(void) const { return AllowedVMStress; } + array GetStressPenaltyParam(void) const { return StressPenaltyParam; } /*! * \brief Get whether a predictor is used for FSI applications. @@ -8608,7 +8577,7 @@ class CConfig { * \brief Get the value of the criteria for applying incremental loading. * \return Value of the log10 of the residual. */ - su2double GetIncLoad_Criteria(unsigned short val_var) const { return IncLoad_Criteria[val_var]; } + su2double GetIncLoad_Criteria(unsigned short val_var) const { return inc_crit[val_var]; } /*! * \brief Get the relaxation method chosen for the simulation @@ -9062,13 +9031,13 @@ class CConfig { * \brief Get the length of the analytic RECTANGLE or BOX grid in the specified coordinate direction. * \return Length the analytic RECTANGLE or BOX grid in the specified coordinate direction. */ - su2double GetMeshBoxLength(unsigned short val_iDim) const { return Mesh_Box_Length[val_iDim]; } + su2double GetMeshBoxLength(unsigned short val_iDim) const { return mesh_box_length[val_iDim]; } /*! * \brief Get the offset from 0.0 of the analytic RECTANGLE or BOX grid in the specified coordinate direction. * \return Offset from 0.0 the analytic RECTANGLE or BOX grid in the specified coordinate direction. */ - su2double GetMeshBoxOffset(unsigned short val_iDim) const { return Mesh_Box_Offset[val_iDim]; } + su2double GetMeshBoxOffset(unsigned short val_iDim) const { return mesh_box_offset[val_iDim]; } /*! * \brief Get the number of screen output variables requested (maximum 6) diff --git a/Common/include/option_structure.inl b/Common/include/option_structure.inl index ff56e139c01b..e69c7762b27b 100644 --- a/Common/include/option_structure.inl +++ b/Common/include/option_structure.inl @@ -335,25 +335,19 @@ public: }; class COptionDoubleArray : public COptionBase { - su2double * & field; // Reference to the feildname - string name; // identifier for the option - const int size; - su2double * def; - su2double * vals; - su2double * default_value; + string name; // Identifier for the option + const int size; // Number of elements + su2double* field; // Reference to the fieldname public: - COptionDoubleArray(string option_field_name, const int list_size, su2double * & option_field, su2double * default_value) : field(option_field), size(list_size) { - this->name = option_field_name; - this->default_value = default_value; - def = nullptr; - vals = nullptr; + COptionDoubleArray(string option_field_name, const int list_size, su2double* option_field) : + name(option_field_name), + size(list_size), + field(option_field) { } - ~COptionDoubleArray() override { - delete [] def; - delete [] vals; - }; + ~COptionDoubleArray() override {}; + string SetValue(vector option_value) override { COptionBase::SetValue(option_value); // Check that the size is correct @@ -371,27 +365,16 @@ public: newstring.append(" found"); return newstring; } - vals = new su2double[this->size]; for (int i = 0; i < this->size; i++) { istringstream is(option_value[i]); - su2double val; - if (!(is >> val)) { - delete [] vals; + if (!(is >> field[i])) { return badValue(option_value, "su2double array", this->name); } - vals[i] = val; } - this->field = vals; return ""; } - void SetDefault() override { - def = new su2double [size]; - for (int i = 0; i < size; i++) { - def[i] = default_value[i]; - } - this->field = def; - } + void SetDefault() override {} }; class COptionDoubleList : public COptionBase { diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 6d82d0091d8f..405b6f6e26b6 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -365,10 +365,10 @@ void CConfig::addEnumListOption(const string name, unsigned short & input_size, option_map.insert( pair(name, val) ); } -void CConfig::addDoubleArrayOption(const string name, const int size, su2double * & option_field, su2double * default_value) { +void CConfig::addDoubleArrayOption(const string name, const int size, su2double* option_field) { assert(option_map.find(name) == option_map.end()); all_options.insert(pair(name, true)); - COptionBase* val = new COptionDoubleArray(name, size, option_field, default_value); + COptionBase* val = new COptionDoubleArray(name, size, option_field); option_map.insert(pair(name, val)); } @@ -907,24 +907,14 @@ void CConfig::SetPointersNull(void) { Aeroelastic_plunge = nullptr; Aeroelastic_pitch = nullptr; - Velocity_FreeStream = nullptr; - Inc_Velocity_Init = nullptr; CFL_AdaptParam = nullptr; CFL = nullptr; - HTP_Axis = nullptr; PlaneTag = nullptr; - Kappa_Flow = nullptr; - Kappa_AdjFlow = nullptr; - Kappa_Heat = nullptr; - Stations_Bounds = nullptr; ParamDV = nullptr; DV_Value = nullptr; Design_Variable = nullptr; - Hold_GridFixed_Coord = nullptr; - SubsonicEngine_Cyl = nullptr; - EA_IntLimit = nullptr; TimeDOFsADER_DG = nullptr; TimeIntegrationADER_DG = nullptr; WeightsIntegrationADER_DG = nullptr; @@ -946,14 +936,6 @@ void CConfig::SetPointersNull(void) { nKind_SurfaceMovement = 0; Kind_SurfaceMovement = nullptr; LocationStations = nullptr; - Motion_Origin = nullptr; - Translation_Rate = nullptr; - Rotation_Rate = nullptr; - Pitching_Omega = nullptr; - Pitching_Ampl = nullptr; - Pitching_Phase = nullptr; - Plunging_Omega = nullptr; - Plunging_Ampl = nullptr; MarkerMotion_Origin = nullptr; MarkerTranslation_Rate = nullptr; MarkerRotation_Rate = nullptr; @@ -993,12 +975,7 @@ void CConfig::SetPointersNull(void) { RelaxFactorAverage = nullptr; RelaxFactorFourier = nullptr; nSpan_iZones = nullptr; - ExtraRelFacGiles = nullptr; - Mixedout_Coeff = nullptr; - RampRotatingFrame_Coeff = nullptr; - RampOutletPressure_Coeff = nullptr; Kind_TurboMachinery = nullptr; - SineLoad_Coeff = nullptr; Marker_MixingPlaneInterface = nullptr; Marker_TurboBoundIn = nullptr; @@ -1117,9 +1094,9 @@ void CConfig::SetConfig_Options() { addBoolOption("GRAVITY_FORCE", GravityForce, false); /* DESCRIPTION: Apply a body force as a source term (NO, YES) */ addBoolOption("BODY_FORCE", Body_Force, false); - default_body_force[0] = 0.0; default_body_force[1] = 0.0; default_body_force[2] = 0.0; + body_force[0] = 0.0; body_force[1] = 0.0; body_force[2] = 0.0; /* DESCRIPTION: Vector of body force values (BodyForce_X, BodyForce_Y, BodyForce_Z) */ - addDoubleArrayOption("BODY_FORCE_VECTOR", 3, Body_Force_Vector, default_body_force); + addDoubleArrayOption("BODY_FORCE_VECTOR", 3, body_force); /*!\brief RESTART_SOL \n DESCRIPTION: Restart solution from native solution file \n Options: NO, YES \ingroup Config */ addBoolOption("RESTART_SOL", Restart, false); /*!\brief BINARY_RESTART \n DESCRIPTION: Read binary SU2 native restart files. \n Options: YES, NO \ingroup Config */ @@ -1213,11 +1190,11 @@ void CConfig::SetConfig_Options() { /*--- Options related to temperature polynomial coefficients for fluid models. ---*/ /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ - addDoubleArrayOption("CP_POLYCOEFFS", N_POLY_COEFFS, CpPolyCoefficients, default_cp_polycoeffs.data()); + addDoubleArrayOption("CP_POLYCOEFFS", N_POLY_COEFFS, cp_polycoeffs.data()); /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ - addDoubleArrayOption("MU_POLYCOEFFS", N_POLY_COEFFS, MuPolyCoefficients, default_mu_polycoeffs.data()); + addDoubleArrayOption("MU_POLYCOEFFS", N_POLY_COEFFS, mu_polycoeffs.data()); /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ - addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, KtPolyCoefficients, default_kt_polycoeffs.data()); + addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, kt_polycoeffs.data()); /*!\brief REYNOLDS_NUMBER \n DESCRIPTION: Reynolds number (non-dimensional, based on the free-stream values). Needed for viscous solvers. For incompressible solvers the Reynolds length will always be 1.0 \n DEFAULT: 0.0 \ingroup Config */ addDoubleOption("REYNOLDS_NUMBER", Reynolds, 0.0); @@ -1262,8 +1239,8 @@ void CConfig::SetConfig_Options() { /*!\brief INC_DENSITY_INIT \n DESCRIPTION: Initial density for incompressible flows (1.2886 kg/m^3 by default) \ingroup Config*/ addDoubleOption("INC_DENSITY_INIT", Inc_Density_Init, 1.2886); /*!\brief INC_VELOCITY_INIT \n DESCRIPTION: Initial velocity for incompressible flows (1.0,0,0 m/s by default) \ingroup Config*/ - default_vel_inf[0] = 1.0; default_vel_inf[1] = 0.0; default_vel_inf[2] = 0.0; - addDoubleArrayOption("INC_VELOCITY_INIT", 3, Inc_Velocity_Init, default_vel_inf); + vel_init[0] = 1.0; vel_init[1] = 0.0; vel_init[2] = 0.0; + addDoubleArrayOption("INC_VELOCITY_INIT", 3, vel_init); /*!\brief INC_TEMPERATURE_INIT \n DESCRIPTION: Initial temperature for incompressible flows with the energy equation (288.15 K by default) \ingroup Config*/ addDoubleOption("INC_TEMPERATURE_INIT", Inc_Temperature_Init, 288.15); /*!\brief INC_NONDIM \n DESCRIPTION: Non-dimensionalization scheme for incompressible flows. \ingroup Config*/ @@ -1275,9 +1252,9 @@ void CConfig::SetConfig_Options() { /*!\brief INC_OUTLET_DAMPING \n DESCRIPTION: Damping factor applied to the iterative updates to the pressure at a mass flow outlet in incompressible flow (0.1 by default). \ingroup Config*/ addDoubleOption("INC_OUTLET_DAMPING", Inc_Outlet_Damping, 0.1); - default_vel_inf[0] = 1.0; default_vel_inf[1] = 0.0; default_vel_inf[2] = 0.0; + vel_inf[0] = 1.0; vel_inf[1] = 0.0; vel_inf[2] = 0.0; /*!\brief FREESTREAM_VELOCITY\n DESCRIPTION: Free-stream velocity (m/s) */ - addDoubleArrayOption("FREESTREAM_VELOCITY", 3, Velocity_FreeStream, default_vel_inf); + addDoubleArrayOption("FREESTREAM_VELOCITY", 3, vel_inf); /* DESCRIPTION: Free-stream viscosity (1.853E-5 Ns/m^2 (air), 0.798E-3 Ns/m^2 (water)) */ addDoubleOption("FREESTREAM_VISCOSITY", Viscosity_FreeStream, -1.0); /* DESCRIPTION: Thermal conductivity used for heat equation */ @@ -1362,8 +1339,8 @@ void CConfig::SetConfig_Options() { /*--- Options related to various boundary markers ---*/ /*!\brief HTP_AXIS\n DESCRIPTION: Location of the HTP axis*/ - default_htp_axis[0] = 0.0; default_htp_axis[1] = 0.0; - addDoubleArrayOption("HTP_AXIS", 2, HTP_Axis, default_htp_axis); + htp_axis[0] = 0.0; htp_axis[1] = 0.0; + addDoubleArrayOption("HTP_AXIS", 2, htp_axis); /*!\brief MARKER_PLOTTING\n DESCRIPTION: Marker(s) of the surface in the surface flow solution file \ingroup Config*/ addStringListOption("MARKER_PLOTTING", nMarker_Plotting, Marker_Plotting); /*!\brief MARKER_MONITORING\n DESCRIPTION: Marker(s) of the surface where evaluate the non-dimensional coefficients \ingroup Config*/ @@ -1458,8 +1435,8 @@ void CConfig::SetConfig_Options() { addBoolOption("SPATIAL_FOURIER", SpatialFourier, false); /*!\brief GILES_EXTRA_RELAXFACTOR \n DESCRIPTION: the 1st coeff the value of the under relaxation factor to apply to the shroud and hub, * the 2nd coefficient is the the percentage of span-wise height influenced by this extra under relaxation factor.*/ - default_extrarelfac[0] = 0.1; default_extrarelfac[1] = 0.1; - addDoubleArrayOption("GILES_EXTRA_RELAXFACTOR", 2, ExtraRelFacGiles, default_extrarelfac); + extrarelfac[0] = 0.1; extrarelfac[1] = 0.1; + addDoubleArrayOption("GILES_EXTRA_RELAXFACTOR", 2, extrarelfac); /*!\brief AVERAGE_PROCESS_TYPE \n DESCRIPTION: types of mixing process for averaging quantities at the boundaries. \n OPTIONS: see \link MixingProcess_Map \endlink \n DEFAULT: AREA_AVERAGE \ingroup Config*/ addEnumOption("MIXINGPLANE_INTERFACE_KIND", Kind_MixingPlaneInterface, MixingPlaneInterface_Map, NEAREST_SPAN); @@ -1469,25 +1446,25 @@ void CConfig::SetConfig_Options() { /*!\brief PERFORMANCE_AVERAGE_PROCESS_KIND \n DESCRIPTION: types of mixing process for averaging quantities at the boundaries for performance computation. \n OPTIONS: see \link MixingProcess_Map \endlink \n DEFAULT: AREA_AVERAGE \ingroup Config*/ addEnumOption("PERFORMANCE_AVERAGE_PROCESS_KIND", Kind_PerformanceAverageProcess, AverageProcess_Map, AREA); - default_mixedout_coeff[0] = 1.0; default_mixedout_coeff[1] = 1.0E-05; default_mixedout_coeff[2] = 15.0; + mixedout_coeff[0] = 1.0; mixedout_coeff[1] = 1.0E-05; mixedout_coeff[2] = 15.0; /*!\brief MIXEDOUT_COEFF \n DESCRIPTION: the 1st coeff is an under relaxation factor for the Newton method, * the 2nd coefficient is the tolerance for the Newton method, 3rd coefficient is the maximum number of * iteration for the Newton Method.*/ - addDoubleArrayOption("MIXEDOUT_COEFF", 3, Mixedout_Coeff, default_mixedout_coeff); + addDoubleArrayOption("MIXEDOUT_COEFF", 3, mixedout_coeff); /*!\brief RAMP_ROTATING_FRAME\n DESCRIPTION: option to ramp up or down the rotating frame velocity value*/ addBoolOption("RAMP_ROTATING_FRAME", RampRotatingFrame, false); - default_rampRotFrame_coeff[0] = 0; default_rampRotFrame_coeff[1] = 1.0; default_rampRotFrame_coeff[2] = 1000.0; + rampRotFrame_coeff[0] = 0; rampRotFrame_coeff[1] = 1.0; rampRotFrame_coeff[2] = 1000.0; /*!\brief RAMP_ROTATING_FRAME_COEFF \n DESCRIPTION: the 1st coeff is the staring velocity, * the 2nd coeff is the number of iterations for the update, 3rd is the number of iteration */ - addDoubleArrayOption("RAMP_ROTATING_FRAME_COEFF", 3, RampRotatingFrame_Coeff, default_rampRotFrame_coeff); + addDoubleArrayOption("RAMP_ROTATING_FRAME_COEFF", 3, rampRotFrame_coeff); /* DESCRIPTION: AVERAGE_MACH_LIMIT is a limit value for average procedure based on the mass flux. */ addDoubleOption("AVERAGE_MACH_LIMIT", AverageMachLimit, 0.03); /*!\brief RAMP_OUTLET_PRESSURE\n DESCRIPTION: option to ramp up or down the rotating frame velocity value*/ addBoolOption("RAMP_OUTLET_PRESSURE", RampOutletPressure, false); - default_rampOutPres_coeff[0] = 100000.0; default_rampOutPres_coeff[1] = 1.0; default_rampOutPres_coeff[2] = 1000.0; + rampOutPres_coeff[0] = 100000.0; rampOutPres_coeff[1] = 1.0; rampOutPres_coeff[2] = 1000.0; /*!\brief RAMP_OUTLET_PRESSURE_COEFF \n DESCRIPTION: the 1st coeff is the staring outlet pressure, * the 2nd coeff is the number of iterations for the update, 3rd is the number of total iteration till reaching the final outlet pressure value */ - addDoubleArrayOption("RAMP_OUTLET_PRESSURE_COEFF", 3, RampOutletPressure_Coeff, default_rampOutPres_coeff); + addDoubleArrayOption("RAMP_OUTLET_PRESSURE_COEFF", 3, rampOutPres_coeff); /*!\brief MARKER_MIXINGPLANE \n DESCRIPTION: Identify the boundaries in which the mixing plane is applied. \ingroup Config*/ addStringListOption("MARKER_MIXINGPLANE_INTERFACE", nMarker_MixingPlaneInterface, Marker_MixingPlaneInterface); /*!\brief TURBULENT_MIXINGPLANE \n DESCRIPTION: Activate mixing plane also for turbulent quantities \ingroup Config*/ @@ -1545,16 +1522,15 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Actuator disk double surface */ addBoolOption("ACTDISK_SU2_DEF", ActDisk_SU2_DEF, false); /* DESCRIPTION: Definition of the distortion rack (radial number of proves / circumferential density (degree) */ - default_distortion[0] = 5.0; default_distortion[1] = 15.0; - addDoubleArrayOption("DISTORTION_RACK", 2, DistortionRack, default_distortion); + distortion[0] = 5.0; distortion[1] = 15.0; + addDoubleArrayOption("DISTORTION_RACK", 2, distortion); /* DESCRIPTION: Values of the box to impose a subsonic nacellle (mach, Pressure, Temperature) */ - default_eng_val[0]=0.0; default_eng_val[1]=0.0; default_eng_val[2]=0.0; - default_eng_val[3]=0.0; default_eng_val[4]=0.0; - addDoubleArrayOption("SUBSONIC_ENGINE_VALUES", 5, SubsonicEngine_Values, default_eng_val); + eng_val[0]=0.0; eng_val[1]=0.0; eng_val[2]=0.0; eng_val[3]=0.0; eng_val[4]=0.0; + addDoubleArrayOption("SUBSONIC_ENGINE_VALUES", 5, eng_val); /* DESCRIPTION: Coordinates of the box to impose a subsonic nacellle cylinder (Xmin, Ymin, Zmin, Xmax, Ymax, Zmax, Radius) */ - default_eng_cyl[0] = 0.0; default_eng_cyl[1] = 0.0; default_eng_cyl[2] = 0.0; - default_eng_cyl[3] = 1E15; default_eng_cyl[4] = 1E15; default_eng_cyl[5] = 1E15; default_eng_cyl[6] = 1E15; - addDoubleArrayOption("SUBSONIC_ENGINE_CYL", 7, SubsonicEngine_Cyl, default_eng_cyl); + eng_cyl[0] = 0.0; eng_cyl[1] = 0.0; eng_cyl[2] = 0.0; + eng_cyl[3] = 1E15; eng_cyl[4] = 1E15; eng_cyl[5] = 1E15; eng_cyl[6] = 1E15; + addDoubleArrayOption("SUBSONIC_ENGINE_CYL", 7, eng_cyl); /* DESCRIPTION: Engine exhaust boundary marker(s) Format: (nacelle exhaust marker, total nozzle temp, total nozzle pressure, ... )*/ addExhaustOption("MARKER_ENGINE_EXHAUST", nMarker_EngineExhaust, Marker_EngineExhaust, Exhaust_Temperature_Target, Exhaust_Pressure_Target); @@ -1577,9 +1553,9 @@ void CConfig::SetConfig_Options() { addInletOption("MARKER_SINE_LOAD", nMarker_Load_Sine, Marker_Load_Sine, Load_Sine_Amplitude, Load_Sine_Frequency, Load_Sine_Dir); /*!\brief SINE_LOAD\n DESCRIPTION: option to apply the load as a sine*/ addBoolOption("SINE_LOAD", Sine_Load, false); - default_sineload_coeff[0] = 0.0; default_sineload_coeff[1] = 0.0; default_sineload_coeff[2] = 0.0; + sineload_coeff[0] = 0.0; sineload_coeff[1] = 0.0; sineload_coeff[2] = 0.0; /*!\brief SINE_LOAD_COEFF \n DESCRIPTION: the 1st coeff is the amplitude, the 2nd is the frequency, 3rd is the phase in radians */ - addDoubleArrayOption("SINE_LOAD_COEFF", 3, SineLoad_Coeff, default_sineload_coeff); + addDoubleArrayOption("SINE_LOAD_COEFF", 3, sineload_coeff); /*!\brief RAMP_AND_RELEASE\n DESCRIPTION: release the load after applying the ramp*/ addBoolOption("RAMP_AND_RELEASE_LOAD", RampAndRelease, false); @@ -1805,14 +1781,14 @@ void CConfig::SetConfig_Options() { /*!\brief SLOPE_LIMITER_FLOW * DESCRIPTION: Slope limiter for the direct solution. \n OPTIONS: See \link Limiter_Map \endlink \n DEFAULT VENKATAKRISHNAN \ingroup Config*/ addEnumOption("SLOPE_LIMITER_FLOW", Kind_SlopeLimit_Flow, Limiter_Map, VENKATAKRISHNAN); - default_jst_coeff[0] = 0.5; default_jst_coeff[1] = 0.02; + jst_coeff[0] = 0.5; jst_coeff[1] = 0.02; /*!\brief JST_SENSOR_COEFF \n DESCRIPTION: 2nd and 4th order artificial dissipation coefficients for the JST method \ingroup Config*/ - addDoubleArrayOption("JST_SENSOR_COEFF", 2, Kappa_Flow, default_jst_coeff); + addDoubleArrayOption("JST_SENSOR_COEFF", 2, jst_coeff); /*!\brief LAX_SENSOR_COEFF \n DESCRIPTION: 1st order artificial dissipation coefficients for the Lax-Friedrichs method. \ingroup Config*/ addDoubleOption("LAX_SENSOR_COEFF", Kappa_1st_Flow, 0.15); - default_ad_coeff_heat[0] = 0.5; default_ad_coeff_heat[1] = 0.02; + ad_coeff_heat[0] = 0.5; ad_coeff_heat[1] = 0.02; /*!\brief JST_SENSOR_COEFF_HEAT \n DESCRIPTION: 2nd and 4th order artificial dissipation coefficients for the JST method \ingroup Config*/ - addDoubleArrayOption("JST_SENSOR_COEFF_HEAT", 2, Kappa_Heat, default_ad_coeff_heat); + addDoubleArrayOption("JST_SENSOR_COEFF_HEAT", 2, ad_coeff_heat); /*!\brief USE_ACCURATE_FLUX_JACOBIANS \n DESCRIPTION: Use numerically computed Jacobians for AUSM+up(2) and SLAU(2) \ingroup Config*/ addBoolOption("USE_ACCURATE_FLUX_JACOBIANS", Use_Accurate_Jacobians, false); /*!\brief CENTRAL_JACOBIAN_FIX_FACTOR \n DESCRIPTION: Improve the numerical properties (diagonal dominance) of the global Jacobian matrix, 3 to 4 is "optimum" (central schemes) \ingroup Config*/ @@ -1829,9 +1805,9 @@ void CConfig::SetConfig_Options() { /*!\brief SLOPE_LIMITER_ADJFLOW * DESCRIPTION: Slope limiter for the adjoint solution. \n OPTIONS: See \link Limiter_Map \endlink \n DEFAULT VENKATAKRISHNAN \ingroup Config*/ addEnumOption("SLOPE_LIMITER_ADJFLOW", Kind_SlopeLimit_AdjFlow, Limiter_Map, VENKATAKRISHNAN); - default_jst_adj_coeff[0] = 0.5; default_jst_adj_coeff[1] = 0.02; + jst_adj_coeff[0] = 0.5; jst_adj_coeff[1] = 0.02; /*!\brief ADJ_JST_SENSOR_COEFF \n DESCRIPTION: 2nd and 4th order artificial dissipation coefficients for the adjoint JST method. \ingroup Config*/ - addDoubleArrayOption("ADJ_JST_SENSOR_COEFF", 2, Kappa_AdjFlow, default_jst_adj_coeff); + addDoubleArrayOption("ADJ_JST_SENSOR_COEFF", 2, jst_adj_coeff); /*!\brief LAX_SENSOR_COEFF \n DESCRIPTION: 1st order artificial dissipation coefficients for the adjoint Lax-Friedrichs method. \ingroup Config*/ addDoubleOption("ADJ_LAX_SENSOR_COEFF", Kappa_1st_AdjFlow, 0.15); @@ -1884,17 +1860,16 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: parameter for the definition of a complex objective function */ addDoubleOption("DCD_DCMY_VALUE", dCD_dCMy, 0.0); - default_obj_coeff[0]=0.0; default_obj_coeff[1]=0.0; default_obj_coeff[2]=0.0; - default_obj_coeff[3]=0.0; default_obj_coeff[4]=0.0; + obj_coeff[0]=0.0; obj_coeff[1]=0.0; obj_coeff[2]=0.0; obj_coeff[3]=0.0; obj_coeff[4]=0.0; /*!\brief OBJ_CHAIN_RULE_COEFF * \n DESCRIPTION: Coefficients defining the objective function gradient using the chain rule * with area-averaged outlet primitive variables. This is used with the genereralized outflow * objective. \ingroup Config */ - addDoubleArrayOption("OBJ_CHAIN_RULE_COEFF", 5, Obj_ChainRuleCoeff, default_obj_coeff); + addDoubleArrayOption("OBJ_CHAIN_RULE_COEFF", 5, obj_coeff); - default_geo_loc[0] = 0.0; default_geo_loc[1] = 1.0; + geo_loc[0] = 0.0; geo_loc[1] = 1.0; /* DESCRIPTION: Definition of the airfoil section */ - addDoubleArrayOption("GEO_BOUNDS", 2, Stations_Bounds, default_geo_loc); + addDoubleArrayOption("GEO_BOUNDS", 2, geo_loc); /* DESCRIPTION: Identify the body to slice */ addEnumOption("GEO_DESCRIPTION", Geo_Description, Geo_Description_Map, WING); /* DESCRIPTION: Z location of the waterline */ @@ -1903,10 +1878,10 @@ void CConfig::SetConfig_Options() { addUnsignedShortOption("GEO_NUMBER_STATIONS", nWingStations, 25); /* DESCRIPTION: Definition of the airfoil sections */ addDoubleListOption("GEO_LOCATION_STATIONS", nLocationStations, LocationStations); - default_nacelle_location[0] = 0.0; default_nacelle_location[1] = 0.0; default_nacelle_location[2] = 0.0; - default_nacelle_location[3] = 0.0; default_nacelle_location[4] = 0.0; + nacelle_location[0] = 0.0; nacelle_location[1] = 0.0; nacelle_location[2] = 0.0; + nacelle_location[3] = 0.0; nacelle_location[4] = 0.0; /* DESCRIPTION: Definition of the nacelle location (higlite coordinates, tilt angle, toe angle) */ - addDoubleArrayOption("GEO_NACELLE_LOCATION", 5, NacelleLocation, default_nacelle_location); + addDoubleArrayOption("GEO_NACELLE_LOCATION", 5, nacelle_location); /* DESCRIPTION: Output sectional forces for specified markers. */ addBoolOption("GEO_PLOT_STATIONS", Plot_Section_Forces, false); /* DESCRIPTION: Mode of the GDC code (analysis, or gradient) */ @@ -1951,12 +1926,12 @@ void CConfig::SetConfig_Options() { addShortListOption("MESH_BOX_SIZE", nMesh_Box_Size, Mesh_Box_Size); /* DESCRIPTION: List of the length of the RECTANGLE or BOX grid in the x,y,z directions. (default: (1.0,1.0,1.0) ). */ - default_mesh_box_length[0] = 1.0; default_mesh_box_length[1] = 1.0; default_mesh_box_length[2] = 1.0; - addDoubleArrayOption("MESH_BOX_LENGTH", 3, Mesh_Box_Length, default_mesh_box_length); + mesh_box_length[0] = 1.0; mesh_box_length[1] = 1.0; mesh_box_length[2] = 1.0; + addDoubleArrayOption("MESH_BOX_LENGTH", 3, mesh_box_length); /* DESCRIPTION: List of the offset from 0.0 of the RECTANGLE or BOX grid in the x,y,z directions. (default: (0.0,0.0,0.0) ). */ - default_mesh_box_offset[0] = 0.0; default_mesh_box_offset[1] = 0.0; default_mesh_box_offset[2] = 0.0; - addDoubleArrayOption("MESH_BOX_OFFSET", 3, Mesh_Box_Offset, default_mesh_box_offset); + mesh_box_offset[0] = 0.0; mesh_box_offset[1] = 0.0; mesh_box_offset[2] = 0.0; + addDoubleArrayOption("MESH_BOX_OFFSET", 3, mesh_box_offset); /* DESCRIPTION: Determine if the mesh file supports multizone. \n DEFAULT: true (temporarily) */ addBoolOption("MULTIZONE_MESH", Multizone_Mesh, true); @@ -2022,23 +1997,22 @@ void CConfig::SetConfig_Options() { addStringListOption("MARKER_MOVING", nMarker_Moving, Marker_Moving); /* DESCRIPTION: Mach number (non-dimensional, based on the mesh velocity and freestream vals.) */ addDoubleOption("MACH_MOTION", Mach_Motion, 0.0); - default_vel_inf[0] = 0.0; default_vel_inf[1] = 0.0; default_vel_inf[2] = 0.0; /* DESCRIPTION: Coordinates of the rigid motion origin */ - addDoubleArrayOption("MOTION_ORIGIN", 3, Motion_Origin, default_vel_inf); + addDoubleArrayOption("MOTION_ORIGIN", 3, Motion_Origin); /* DESCRIPTION: Translational velocity vector (m/s) in the x, y, & z directions (RIGID_MOTION only) */ - addDoubleArrayOption("TRANSLATION_RATE", 3, Translation_Rate, default_vel_inf); + addDoubleArrayOption("TRANSLATION_RATE", 3, Translation_Rate); /* DESCRIPTION: Angular velocity vector (rad/s) about x, y, & z axes (RIGID_MOTION only) */ - addDoubleArrayOption("ROTATION_RATE", 3, Rotation_Rate, default_vel_inf); + addDoubleArrayOption("ROTATION_RATE", 3, Rotation_Rate); /* DESCRIPTION: Pitching angular freq. (rad/s) about x, y, & z axes (RIGID_MOTION only) */ - addDoubleArrayOption("PITCHING_OMEGA", 3, Pitching_Omega, default_vel_inf); + addDoubleArrayOption("PITCHING_OMEGA", 3, Pitching_Omega); /* DESCRIPTION: Pitching amplitude (degrees) about x, y, & z axes (RIGID_MOTION only) */ - addDoubleArrayOption("PITCHING_AMPL", 3, Pitching_Ampl, default_vel_inf); + addDoubleArrayOption("PITCHING_AMPL", 3, Pitching_Ampl); /* DESCRIPTION: Pitching phase offset (degrees) about x, y, & z axes (RIGID_MOTION only) */ - addDoubleArrayOption("PITCHING_PHASE", 3, Pitching_Phase, default_vel_inf); + addDoubleArrayOption("PITCHING_PHASE", 3, Pitching_Phase); /* DESCRIPTION: Plunging angular freq. (rad/s) in x, y, & z directions (RIGID_MOTION only) */ - addDoubleArrayOption("PLUNGING_OMEGA", 3, Plunging_Omega, default_vel_inf); + addDoubleArrayOption("PLUNGING_OMEGA", 3, Plunging_Omega); /* DESCRIPTION: Plunging amplitude (m) in x, y, & z directions (RIGID_MOTION only) */ - addDoubleArrayOption("PLUNGING_AMPL", 3, Plunging_Ampl, default_vel_inf); + addDoubleArrayOption("PLUNGING_AMPL", 3, Plunging_Ampl); /* DESCRIPTION: Coordinates of the rigid motion origin */ addDoubleListOption("SURFACE_MOTION_ORIGIN", nMarkerMotion_Origin, MarkerMotion_Origin); /* DESCRIPTION: Translational velocity vector (m/s) in the x, y, & z directions (DEFORMING only) */ @@ -2128,9 +2102,9 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Evaluate equivalent area on the Near-Field */ addBoolOption("EQUIV_AREA", EquivArea, false); - default_ea_lim[0] = 0.0; default_ea_lim[1] = 1.0; default_ea_lim[2] = 1.0; + ea_lim[0] = 0.0; ea_lim[1] = 1.0; ea_lim[2] = 1.0; /* DESCRIPTION: Integration limits of the equivalent area ( xmin, xmax, Dist_NearField ) */ - addDoubleArrayOption("EA_INT_LIMIT", 3, EA_IntLimit, default_ea_lim); + addDoubleArrayOption("EA_INT_LIMIT", 3, ea_lim); /* DESCRIPTION: Equivalent area scaling factor */ addDoubleOption("EA_SCALE_FACTOR", EA_ScaleFactor, 1.0); @@ -2176,10 +2150,10 @@ void CConfig::SetConfig_Options() { addEnumOption("DV_SENSITIVITY_FORMAT", Sensitivity_FileFormat, Sensitivity_Map, SU2_NATIVE); /* DESCRIPTION: Hold the grid fixed in a region */ addBoolOption("HOLD_GRID_FIXED", Hold_GridFixed, false); - default_grid_fix[0] = -1E15; default_grid_fix[1] = -1E15; default_grid_fix[2] = -1E15; - default_grid_fix[3] = 1E15; default_grid_fix[4] = 1E15; default_grid_fix[5] = 1E15; + grid_fix[0] = -1E15; grid_fix[1] = -1E15; grid_fix[2] = -1E15; + grid_fix[3] = 1E15; grid_fix[4] = 1E15; grid_fix[5] = 1E15; /* DESCRIPTION: Coordinates of the box where the grid will be deformed (Xmin, Ymin, Zmin, Xmax, Ymax, Zmax) */ - addDoubleArrayOption("HOLD_GRID_FIXED_COORD", 6, Hold_GridFixed_Coord, default_grid_fix); + addDoubleArrayOption("HOLD_GRID_FIXED_COORD", 6, grid_fix); /*!\par CONFIG_CATEGORY: Deformable mesh \ingroup Config*/ /*--- option related to deformable meshes ---*/ @@ -2294,8 +2268,8 @@ void CConfig::SetConfig_Options() { /*!\brief REFERENCE_NODE_PENALTY\n DESCRIPTION: Penalty weight value for the objective function \ingroup Config*/ addDoubleOption("REFERENCE_NODE_PENALTY", RefNode_Penalty, 1E3); - /*!\brief ALLOWED_VONMISSES_STRESS\n DESCRIPTION: Maximum allowed stress for structural optimization \ingroup Config*/ - addDoubleOption("ALLOWED_VONMISSES_STRESS", AllowedVMStress, 1.0); + /*!\brief STRESS_PENALTY_PARAM\n DESCRIPTION: Maximum allowed stress and KS exponent for structural optimization \ingroup Config*/ + addDoubleArrayOption("STRESS_PENALTY_PARAM", 2, StressPenaltyParam.data()); /*!\brief REGIME_TYPE \n DESCRIPTION: Geometric condition \n OPTIONS: see \link Struct_Map \endlink \ingroup Config*/ addEnumOption("GEOMETRIC_CONDITIONS", Kind_Struct_Solver, Struct_Map, SMALL_DEFORMATIONS); @@ -2345,9 +2319,9 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Maximum number of increments of the */ addUnsignedLongOption("NUMBER_INCREMENTS", IncLoad_Nincrements, 10); - default_inc_crit[0] = 0.0; default_inc_crit[1] = 0.0; default_inc_crit[2] = 0.0; + inc_crit[0] = 0.0; inc_crit[1] = 0.0; inc_crit[2] = 0.0; /* DESCRIPTION: Definition of the UTOL RTOL ETOL*/ - addDoubleArrayOption("INCREMENTAL_CRITERIA", 3, IncLoad_Criteria, default_inc_crit); + addDoubleArrayOption("INCREMENTAL_CRITERIA", 3, inc_crit); /* DESCRIPTION: Use of predictor */ addBoolOption("PREDICTOR", Predictor, false); @@ -2486,11 +2460,11 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Rotation of the volumetric heat source respect to Z axis */ addDoubleOption("HEAT_SOURCE_ROTATION_Z", Heat_Source_Rot_Z, 0.0); /* DESCRIPTION: Position of heat source center (Heat_Source_Center_X, Heat_Source_Center_Y, Heat_Source_Center_Z) */ - default_hs_center[0] = 0.0; default_hs_center[1] = 0.0; default_hs_center[2] = 0.0; - addDoubleArrayOption("HEAT_SOURCE_CENTER", 3, Heat_Source_Center, default_hs_center); + hs_center[0] = 0.0; hs_center[1] = 0.0; hs_center[2] = 0.0; + addDoubleArrayOption("HEAT_SOURCE_CENTER", 3, hs_center); /* DESCRIPTION: Vector of heat source radii (Heat_Source_Axes_A, Heat_Source_Axes_B, Heat_Source_Axes_C) */ - default_hs_axes[0] = 1.0; default_hs_axes[1] = 1.0; default_hs_axes[2] = 1.0; - addDoubleArrayOption("HEAT_SOURCE_AXES", 3, Heat_Source_Axes, default_hs_axes); + hs_axes[0] = 1.0; hs_axes[1] = 1.0; hs_axes[2] = 1.0; + addDoubleArrayOption("HEAT_SOURCE_AXES", 3, hs_axes); /*!\brief MARKER_EMISSIVITY DESCRIPTION: Wall emissivity of the marker for radiation purposes \n * Format: ( marker, emissivity of the marker, ... ) \ingroup Config */ @@ -2555,8 +2529,8 @@ void CConfig::SetConfig_Options() { addEnumOption("FFD_COORD_SYSTEM", FFD_CoordSystem, CoordSystem_Map, CARTESIAN); /* DESCRIPTION: Axis information for the spherical and cylindrical coord system */ - default_ffd_axis[0] = 0.0; default_ffd_axis[1] = 0.0; default_ffd_axis[2] =0.0; - addDoubleArrayOption("FFD_AXIS", 3, FFD_Axis, default_ffd_axis); + ffd_axis[0] = 0.0; ffd_axis[1] = 0.0; ffd_axis[2] =0.0; + addDoubleArrayOption("FFD_AXIS", 3, ffd_axis); /* DESCRIPTION: Number of total iterations in the FFD point inversion */ addUnsignedShortOption("FFD_ITERATIONS", nFFD_Iter, 500); @@ -2595,8 +2569,8 @@ void CConfig::SetConfig_Options() { addEnumOption("FFD_BLENDING", FFD_Blending, Blending_Map, BEZIER ); /* DESCRIPTION: Order of the BSplines for BSpline Blending function */ - default_ffd_coeff[0] = 2; default_ffd_coeff[1] = 2; default_ffd_coeff[2] = 2; - addDoubleArrayOption("FFD_BSPLINE_ORDER", 3, FFD_BSpline_Order, default_ffd_coeff); + ffd_coeff[0] = 2; ffd_coeff[1] = 2; ffd_coeff[2] = 2; + addDoubleArrayOption("FFD_BSPLINE_ORDER", 3, ffd_coeff); /*--- Options for the automatic differentiation methods ---*/ /*!\par CONFIG_CATEGORY: Automatic Differentation options\ingroup Config*/ @@ -3265,18 +3239,18 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ /*--- Compute x-velocity with a safegaurd for 0.0. ---*/ su2double Vx = 1e-10; - if (Inc_Velocity_Init[0] != 0.0) { - Vx = Inc_Velocity_Init[0]; + if (vel_init[0] != 0.0) { + Vx = vel_init[0]; } /*--- Compute the angle-of-attack and sideslip. ---*/ su2double alpha = 0.0, beta = 0.0; if (val_nDim == 2) { - alpha = atan(Inc_Velocity_Init[1]/Vx)*180.0/PI_NUMBER; + alpha = atan(vel_init[1]/Vx)*180.0/PI_NUMBER; } else { - alpha = atan(Inc_Velocity_Init[2]/Vx)*180.0/PI_NUMBER; - beta = atan(Inc_Velocity_Init[1]/Vx)*180.0/PI_NUMBER; + alpha = atan(vel_init[2]/Vx)*180.0/PI_NUMBER; + beta = atan(vel_init[1]/Vx)*180.0/PI_NUMBER; } /*--- Set alpha and beta in the config class. ---*/ @@ -3796,28 +3770,28 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ if(GetGrid_Movement() && RampRotatingFrame && !DiscreteAdjoint){ FinalRotation_Rate_Z = Rotation_Rate[2]; if(abs(FinalRotation_Rate_Z) > 0.0){ - Rotation_Rate[2] = RampRotatingFrame_Coeff[0]; + Rotation_Rate[2] = rampRotFrame_coeff[0]; } } if(RampOutletPressure && !DiscreteAdjoint){ for (iMarker = 0; iMarker < nMarker_Giles; iMarker++){ if (Kind_Data_Giles[iMarker] == STATIC_PRESSURE || Kind_Data_Giles[iMarker] == STATIC_PRESSURE_1D || Kind_Data_Giles[iMarker] == RADIAL_EQUILIBRIUM ){ - FinalOutletPressure = Giles_Var1[iMarker]; - Giles_Var1[iMarker] = RampOutletPressure_Coeff[0]; + FinalOutletPressure = Giles_Var1[iMarker]; + Giles_Var1[iMarker] = rampOutPres_coeff[0]; } } for (iMarker = 0; iMarker < nMarker_Riemann; iMarker++){ if (Kind_Data_Riemann[iMarker] == STATIC_PRESSURE || Kind_Data_Riemann[iMarker] == RADIAL_EQUILIBRIUM){ - FinalOutletPressure = Riemann_Var1[iMarker]; - Riemann_Var1[iMarker] = RampOutletPressure_Coeff[0]; - } + FinalOutletPressure = Riemann_Var1[iMarker]; + Riemann_Var1[iMarker] = rampOutPres_coeff[0]; } } + } /*--- Check on extra Relaxation factor for Giles---*/ - if(ExtraRelFacGiles[1] > 0.5){ - ExtraRelFacGiles[1] = 0.5; + if(extrarelfac[1] > 0.5){ + extrarelfac[1] = 0.5; } /*--- Use the various rigid-motion input frequencies to determine the period to be used with harmonic balance cases. There are THREE types of motion to consider, namely: rotation, pitching, and plunging. @@ -4001,12 +3975,12 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ Kind_Solver == FEM_EULER) Kind_Turb_Model = NONE; - Kappa_2nd_Flow = Kappa_Flow[0]; - Kappa_4th_Flow = Kappa_Flow[1]; - Kappa_2nd_AdjFlow = Kappa_AdjFlow[0]; - Kappa_4th_AdjFlow = Kappa_AdjFlow[1]; - Kappa_2nd_Heat = Kappa_Heat[0]; - Kappa_4th_Heat = Kappa_Heat[1]; + Kappa_2nd_Flow = jst_coeff[0]; + Kappa_4th_Flow = jst_coeff[1]; + Kappa_2nd_AdjFlow = jst_adj_coeff[0]; + Kappa_4th_AdjFlow = jst_adj_coeff[1]; + Kappa_2nd_Heat = ad_coeff_heat[0]; + Kappa_4th_Heat = ad_coeff_heat[1]; /*--- Make the MG_PreSmooth, MG_PostSmooth, and MG_CorrecSmooth arrays consistent with nMGLevels ---*/ @@ -4173,8 +4147,8 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ Kind_ConvNumScheme_Flow = Kind_ConvNumScheme_AdjFlow; Kind_Centered_Flow = Kind_Centered_AdjFlow; Kind_Upwind_Flow = Kind_Upwind_AdjFlow; - Kappa_Flow[0] = Kappa_AdjFlow[0]; - Kappa_Flow[1] = Kappa_AdjFlow[1]; + Kappa_2nd_Flow = jst_adj_coeff[0]; + Kappa_4th_Flow = jst_adj_coeff[1]; } if (Update_AoA_Iter_Limit == 0 && Fixed_CL_Mode) { @@ -4332,8 +4306,8 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ for (unsigned short iSections = 0; iSections < nLocationStations; iSections++) { LocationStations[iSections] += EPS; } - Stations_Bounds[0] += EPS; - Stations_Bounds[1] += EPS; + geo_loc[0] += EPS; + geo_loc[1] += EPS; } /*--- Length based parameter for slope limiters uses a default value of @@ -4367,26 +4341,19 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ Highlite_Area = Highlite_Area/144.0; SemiSpan = SemiSpan/12.0; - EA_IntLimit[0] = EA_IntLimit[0]/12.0; - EA_IntLimit[1] = EA_IntLimit[1]/12.0; - EA_IntLimit[2] = EA_IntLimit[2]/12.0; + ea_lim[0] /= 12.0; + ea_lim[1] /= 12.0; + ea_lim[2] /= 12.0; if (Geo_Description != NACELLE) { for (unsigned short iSections = 0; iSections < nLocationStations; iSections++) { LocationStations[iSections] = LocationStations[iSections]/12.0; } - Stations_Bounds[0] = Stations_Bounds[0]/12.0; - Stations_Bounds[1] = Stations_Bounds[1]/12.0; + geo_loc[0] /= 12.0; + geo_loc[1] /= 12.0; } - SubsonicEngine_Cyl[0] = SubsonicEngine_Cyl[0]/12.0; - SubsonicEngine_Cyl[1] = SubsonicEngine_Cyl[1]/12.0; - SubsonicEngine_Cyl[2] = SubsonicEngine_Cyl[2]/12.0; - SubsonicEngine_Cyl[3] = SubsonicEngine_Cyl[3]/12.0; - SubsonicEngine_Cyl[4] = SubsonicEngine_Cyl[4]/12.0; - SubsonicEngine_Cyl[5] = SubsonicEngine_Cyl[5]/12.0; - SubsonicEngine_Cyl[6] = SubsonicEngine_Cyl[6]/12.0; - + for (int i=0; i<7; ++i) eng_cyl[i] /= 12.0; } if ((Kind_Turb_Model != SA) && (Kind_Trans_Model == BC)){ @@ -5635,14 +5602,14 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) { } if (Fixed_CM_Mode) { cout << "Fixed CM mode, target value: " << Target_CM << "." << endl; - cout << "HTP rotation axis (X,Z): ("<< HTP_Axis[0] <<", "<< HTP_Axis[1] <<")."<< endl; + cout << "HTP rotation axis (X,Z): ("<< htp_axis[0] <<", "<< htp_axis[1] <<")."<< endl; } } if (EquivArea) { cout <<"The equivalent area is going to be evaluated on the near-field."<< endl; - cout <<"The lower integration limit is "<GetKind_GridMovement() == AEROELASTIC_RIGID_MOTION) { su2double Omega, dt, psi; dt = config->GetDelta_UnstTimeND(); - Omega = (config->GetRotation_Rate(3)/config->GetOmega_Ref()); + Omega = config->GetRotation_Rate(2)/config->GetOmega_Ref(); psi = Omega*(dt*TimeIter); /*--- Correct for the airfoil starting position (This is hardcoded in here) ---*/ diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 83221d895ccc..fae1a2d6fda2 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -5841,7 +5841,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, unsigned short iDim, iVar, jVar, kVar; unsigned long iVertex, iPoint, Point_Normal; - su2double P_Total, T_Total, P_static, T_static, Rho_static, *Mach, *Flow_Dir, Area, UnitNormal[3]; + const su2double *Flow_Dir, *Mach; + su2double P_Total, T_Total, P_static, T_static, Rho_static, Area, UnitNormal[MAXNDIM]; su2double *Velocity_b, Velocity2_b, Enthalpy_b, Energy_b, StaticEnergy_b, Density_b, Kappa_b, Chi_b, Pressure_b, Temperature_b; su2double *Velocity_e, Velocity2_e, VelMag_e, Enthalpy_e, Entropy_e, Energy_e = 0.0, StaticEnthalpy_e, StaticEnergy_e, Density_e = 0.0, Pressure_e; su2double *Velocity_i, Velocity2_i, Enthalpy_i, Energy_i, StaticEnergy_i, Density_i, Kappa_i, Chi_i, Pressure_i, SoundSpeed_i; @@ -6352,7 +6353,8 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain unsigned short iDim, iVar, jVar, kVar, iSpan; unsigned long iPoint, Point_Normal, oldVertex, iVertex; - su2double P_Total, T_Total, *Flow_Dir; + const su2double *Flow_Dir; + su2double P_Total, T_Total; su2double *Velocity_b, Velocity2_b, Enthalpy_b, Energy_b, StaticEnergy_b, Density_b, Kappa_b, Chi_b, Pressure_b, Temperature_b; su2double *Velocity_e, Velocity2_e, Enthalpy_e, Entropy_e, Energy_e = 0.0, StaticEnthalpy_e, StaticEnergy_e, Density_e = 0.0, Pressure_e; su2double *Velocity_i, Velocity2_i, Enthalpy_i, Energy_i, StaticEnergy_i, Density_i, Kappa_i, Chi_i, Pressure_i, SoundSpeed_i; @@ -7036,7 +7038,7 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu su2double relfacFouCfg = config->GetGiles_RelaxFactorFourier(Marker_Tag); su2double *Normal; su2double TwoPiThetaFreq_Pitch, pitch,theta; - const su2double *SpanWiseValues = nullptr; + const su2double *SpanWiseValues = nullptr, *FlowDir; su2double spanPercent, extrarelfacAvg = 0.0, deltaSpan = 0.0, relfacAvg, relfacFou, coeffrelfacAvg = 0.0; unsigned short Turbo_Flag; @@ -7053,7 +7055,7 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu S_boundary = new su2double[8]; su2double AvgMach , *cj, GilesBeta, *delta_c, **R_Matrix, *deltaprim, **R_c_inv,**R_c, alphaIn_BC, gammaIn_BC = 0, - P_Total, T_Total, *FlowDir, Enthalpy_BC, Entropy_BC, *R, *c_avg,*dcjs, Beta_inf2, c2js_Re, c3js_Re, cOutjs_Re, avgVel2 =0.0; + P_Total, T_Total, Enthalpy_BC, Entropy_BC, *R, *c_avg,*dcjs, Beta_inf2, c2js_Re, c3js_Re, cOutjs_Re, avgVel2 =0.0; long freq; @@ -8233,7 +8235,7 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con unsigned long iVertex, iPoint; su2double *V_inlet, *V_domain; - su2double Density, Pressure, Temperature, Energy, *Vel, Velocity2; + su2double Density, Energy, Velocity2; su2double Gas_Constant = config->GetGas_ConstantND(); bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); @@ -8246,9 +8248,9 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con so all flow variables can be imposed at the inlet. First, retrieve the specified values for the primitive variables. ---*/ - Temperature = config->GetInlet_Temperature(Marker_Tag); - Pressure = config->GetInlet_Pressure(Marker_Tag); - Vel = config->GetInlet_Velocity(Marker_Tag); + auto Temperature = config->GetInlet_Temperature(Marker_Tag); + auto Pressure = config->GetInlet_Pressure(Marker_Tag); + auto Vel = config->GetInlet_Velocity(Marker_Tag); /*--- Non-dim. the inputs if necessary. ---*/ @@ -10202,7 +10204,7 @@ void CEulerSolver::SetFreeStream_TurboSolution(CConfig *config) { unsigned long iPoint; unsigned short iDim; unsigned short iZone = config->GetiZone(); - su2double *turboVelocity, *cartVelocity, *turboNormal; + su2double *turboVelocity, *cartVelocity; su2double Alpha = config->GetAoA()*PI_NUMBER/180.0; su2double Mach = config->GetMach(); @@ -10211,7 +10213,7 @@ void CEulerSolver::SetFreeStream_TurboSolution(CConfig *config) { turboVelocity = new su2double[nDim]; cartVelocity = new su2double[nDim]; - turboNormal = config->GetFreeStreamTurboNormal(); + auto turboNormal = config->GetFreeStreamTurboNormal(); GetFluidModel()->SetTDState_Prho(Pressure_Inf, Density_Inf); SoundSpeed = GetFluidModel()->GetSoundSpeed(); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index d09225cf6850..c25fcd582ebd 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -1329,9 +1329,12 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, const bool topology_mode = config->GetTopology_Optimization(); const auto simp_exponent = config->GetSIMP_Exponent(); + const auto stressParam = config->GetStressPenaltyParam(); + const su2double stress_scale = 1.0 / stressParam[0]; + const su2double ks_mult = stressParam[1]; + const unsigned short nStress = (nDim == 2) ? 3 : 6; - const su2double stressScale = 1.0 / config->GetAllowedVMStress(); su2double StressPenalty = 0.0; su2double MaxVonMises_Stress = 0.0; @@ -1412,7 +1415,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, auto elStress = numerics[NUM_TERM]->Compute_Averaged_NodalStress(element, config); - stressPen += pow(max(0.0, elStress*simp_penalty*stressScale - 1.0), 2); + stressPen += exp(ks_mult * elStress*simp_penalty*stress_scale); wasActive = AD::BeginPassive(); @@ -1469,6 +1472,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, /*--- Reduce the stress penalty over all ranks ---*/ SU2_MPI::Allreduce(&StressPenalty, &Total_OFStressPenalty, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + Total_OFStressPenalty = log(Total_OFStressPenalty)/ks_mult - 1.0; bool outputReactions = false; diff --git a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp index d923e4fadf89..df70ceb8dc43 100644 --- a/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp +++ b/SU2_CFD/src/solvers/CFEM_DG_EulerSolver.cpp @@ -7721,7 +7721,7 @@ void CFEM_DG_EulerSolver::BoundaryStates_Riemann(CConfig *confi pressure and temperature as well as the flow direction. */ su2double P_Total = config->GetRiemann_Var1(Marker_Tag); su2double T_Total = config->GetRiemann_Var2(Marker_Tag); - su2double *Flow_Dir = config->GetRiemann_FlowDir(Marker_Tag); + auto Flow_Dir = config->GetRiemann_FlowDir(Marker_Tag); P_Total /= config->GetPressure_Ref(); T_Total /= config->GetTemperature_Ref(); @@ -7782,7 +7782,7 @@ void CFEM_DG_EulerSolver::BoundaryStates_Riemann(CConfig *confi temperature as well as the three components of the Mach number. */ su2double P_static = config->GetRiemann_Var1(Marker_Tag); su2double T_static = config->GetRiemann_Var2(Marker_Tag); - su2double *Mach = config->GetRiemann_FlowDir(Marker_Tag); + auto Mach = config->GetRiemann_FlowDir(Marker_Tag); P_static /= config->GetPressure_Ref(); T_static /= config->GetTemperature_Ref(); @@ -7833,7 +7833,7 @@ void CFEM_DG_EulerSolver::BoundaryStates_Riemann(CConfig *confi temperature as well as the three components of the Mach number. */ su2double P_static = config->GetRiemann_Var1(Marker_Tag); su2double Rho_static = config->GetRiemann_Var2(Marker_Tag); - su2double *Mach = config->GetRiemann_FlowDir(Marker_Tag); + auto Mach = config->GetRiemann_FlowDir(Marker_Tag); P_static /= config->GetPressure_Ref(); Rho_static /= config->GetDensity_Ref(); @@ -7883,7 +7883,7 @@ void CFEM_DG_EulerSolver::BoundaryStates_Riemann(CConfig *confi flow direction. Retrieve the non-dimensional data. */ su2double Density_e = config->GetRiemann_Var1(Marker_Tag); su2double VelMag_e = config->GetRiemann_Var2(Marker_Tag); - su2double *Flow_Dir = config->GetRiemann_FlowDir(Marker_Tag); + auto Flow_Dir = config->GetRiemann_FlowDir(Marker_Tag); Density_e /= config->GetDensity_Ref(); VelMag_e /= config->GetVelocity_Ref(); diff --git a/SU2_PY/SU2/io/historyMap.py b/SU2_PY/SU2/io/historyMap.py index 157218eaeaaf..f809b9a426dc 100644 --- a/SU2_PY/SU2/io/historyMap.py +++ b/SU2_PY/SU2/io/historyMap.py @@ -374,6 +374,10 @@ 'GROUP': 'D_ENGINE_OUTPUT', 'HEADER': 'd[SolidCDrag]', 'TYPE': 'D_COEFFICIENT'}, + 'D_STRESS_PENALTY': {'DESCRIPTION': 'Derivative value', + 'GROUP': 'D_STRUCT_COEFF', + 'HEADER': 'd[StressPen]', + 'TYPE': 'D_COEFFICIENT'}, 'D_SURFACE_MACH': {'DESCRIPTION': 'Derivative value', 'GROUP': 'D_FLOW_COEFF', 'HEADER': 'd[Avg_Mach]', @@ -868,6 +872,10 @@ 'GROUP': 'ENGINE_OUTPUT', 'HEADER': 'SolidCDrag', 'TYPE': 'COEFFICIENT'}, + 'STRESS_PENALTY': {'DESCRIPTION': '', + 'GROUP': 'STRUCT_COEFF', + 'HEADER': 'StressPen', + 'TYPE': 'COEFFICIENT'}, 'SURFACE_MACH': {'DESCRIPTION': 'Total average mach number on all markers set ' 'in MARKER_ANALYZE', 'GROUP': 'FLOW_COEFF', @@ -1114,6 +1122,11 @@ 'GROUP': 'TAVG_D_ENGINE_OUTPUT', 'HEADER': 'dtavg[SolidCDrag]', 'TYPE': 'TAVG_D_COEFFICIENT'}, + 'TAVG_D_STRESS_PENALTY': {'DESCRIPTION': 'weighted time average derivative ' + 'value', + 'GROUP': 'TAVG_D_STRUCT_COEFF', + 'HEADER': 'dtavg[StressPen]', + 'TYPE': 'TAVG_D_COEFFICIENT'}, 'TAVG_D_SURFACE_MACH': {'DESCRIPTION': 'weighted time average derivative ' 'value', 'GROUP': 'TAVG_D_FLOW_COEFF', @@ -1293,6 +1306,10 @@ 'GROUP': 'TAVG_ENGINE_OUTPUT', 'HEADER': 'tavg[SolidCDrag]', 'TYPE': 'TAVG_COEFFICIENT'}, + 'TAVG_STRESS_PENALTY': {'DESCRIPTION': 'weighted time average value', + 'GROUP': 'TAVG_STRUCT_COEFF', + 'HEADER': 'tavg[StressPen]', + 'TYPE': 'TAVG_COEFFICIENT'}, 'TAVG_SURFACE_MACH': {'DESCRIPTION': 'weighted time average value', 'GROUP': 'TAVG_FLOW_COEFF', 'HEADER': 'tavg[Avg_Mach]', diff --git a/TestCases/fea_topology/config.cfg b/TestCases/fea_topology/config.cfg index ec6229018ba0..94c1e956f3f8 100644 --- a/TestCases/fea_topology/config.cfg +++ b/TestCases/fea_topology/config.cfg @@ -50,12 +50,16 @@ TOPOL_OPTIM_OUTFILE= grad_ref_node.dat % - VOLUME_FRACTION: Volume average of rho, 1 means totally solid; % - TOPOL_DISCRETENESS: Volume average of 4*rho*(1-rho), % 0 means perfect solid-void topology. +% - STRESS_PENALTY: KS-aggregated maximum element-average VM stress, +% use it as a <= 0 constraint. % REFERENCE_NODE can be used in lieu of compliance for simple load cases. OBJECTIVE_FUNCTION= REFERENCE_NODE REFERENCE_NODE= 5225 REFERENCE_NODE_DISPLACEMENT= (0.0, 0.0) REFERENCE_NODE_PENALTY= 1.0 DESIGN_VARIABLE_FEA= YOUNG_MODULUS +% Parameters for the corresponding OF (allowed stress and KS multiplier). +STRESS_PENALTY_PARAM= (1.0, 10.0) % ITER=1 % From aa3fa3f05013bc6ae62dcab425239edd71dfc833 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 22 Jan 2021 18:29:22 +0000 Subject: [PATCH 4/9] remove more SU2_MSH --- Common/include/CConfig.hpp | 53 +-------------- Common/include/option_structure.hpp | 35 ---------- Common/src/CConfig.cpp | 100 ++-------------------------- SU2_CFD/src/output/CFlowOutput.cpp | 4 +- 4 files changed, 8 insertions(+), 184 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 0ea1ffc0c9e7..3985d25e88ab 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -109,6 +109,7 @@ class CConfig { su2double Opt_RelaxFactor; /*!< \brief Scale factor for the line search. */ su2double Opt_LineSearch_Bound; /*!< \brief Bounds for the line search. */ su2double StartTime; + unsigned short SmoothNumGrid; /*!< \brief Smooth the numerical grid. */ bool ContinuousAdjoint, /*!< \brief Flag to know if the code is solving an adjoint problem. */ Viscous, /*!< \brief Flag to know if the code is solving a viscous problem. */ EquivArea, /*!< \brief Flag to know if the code is going to compute and plot the equivalent area. */ @@ -123,8 +124,6 @@ class CConfig { Low_Mach_Precon, /*!< \brief Flag to know if we are using a low Mach number preconditioner. */ Low_Mach_Corr, /*!< \brief Flag to know if we are using a low Mach number correction. */ GravityForce, /*!< \brief Flag to know if the gravity force is incuded in the formulation. */ - SmoothNumGrid, /*!< \brief Smooth the numerical grid. */ - AdaptBoundary, /*!< \brief Adapt the elements on the boundary. */ SubsonicEngine, /*!< \brief Engine intake subsonic region. */ Frozen_Visc_Cont, /*!< \brief Flag for cont. adjoint problem with/without frozen viscosity. */ Frozen_Visc_Disc, /*!< \brief Flag for disc. adjoint problem with/without frozen viscosity. */ @@ -167,10 +166,8 @@ class CConfig { unsigned short Continuous_Eqns; /*!< \brief Which equations to treat continuously (Hybrid adjoint)*/ unsigned short Discrete_Eqns; /*!< \brief Which equations to treat discretely (Hybrid adjoint). */ unsigned short *Design_Variable; /*!< \brief Kind of design variable. */ - unsigned short Kind_Adaptation; /*!< \brief Kind of numerical grid adaptation. */ unsigned short nTimeInstances; /*!< \brief Number of periodic time instances for harmonic balance. */ su2double HarmonicBalance_Period; /*!< \brief Period of oscillation to be used with harmonic balance computations. */ - su2double New_Elem_Adapt; /*!< \brief Elements to adapt in the numerical grid adaptation process. */ su2double Delta_UnstTime, /*!< \brief Time step for unsteady computations. */ Delta_UnstTimeND; /*!< \brief Time step for unsteady computations (non dimensional). */ su2double Delta_DynTime, /*!< \brief Time step for dynamic structural computations. */ @@ -712,9 +709,7 @@ class CConfig { *Marker_CfgFile_DV, /*!< \brief Global index for design variable markers using the config information. */ *Marker_CfgFile_PerBound; /*!< \brief Global index for periodic boundaries using the config information. */ string *PlaneTag; /*!< \brief Global index for the plane adaptation (upper, lower). */ - su2double DualVol_Power; /*!< \brief Power for the dual volume in the grid adaptation sensor. */ su2double *nBlades; /*!< \brief number of blades for turbomachinery computation. */ - unsigned short Analytical_Surface; /*!< \brief Information about the analytical definition of the surface for grid adaptation. */ unsigned short Geo_Description; /*!< \brief Description of the geometry. */ unsigned short Mesh_FileFormat; /*!< \brief Mesh input format. */ unsigned short Tab_FileFormat; /*!< \brief Format of the output files. */ @@ -1458,20 +1453,6 @@ class CConfig { */ const su2double *GetDistortionRack(void) const { return distortion; } - /*! - * \brief Get the power of the dual volume in the grid adaptation sensor. - * \return Power of the dual volume in the grid adaptation sensor. - */ - su2double GetDualVol_Power(void) const { return DualVol_Power; } - - /*! - * \brief Get Information about if there is an analytical definition of the surface for doing the - * grid adaptation. - * \return Definition of the surfaces. NONE implies that there isn't any analytical definition - * and it will use and interpolation. - */ - unsigned short GetAnalytical_Surface(void) const { return Analytical_Surface; } - /*! * \brief Get Description of the geometry to be analyzed */ @@ -4157,18 +4138,6 @@ class CConfig { */ unsigned short GetKind_SGS_Model(void) const { return Kind_SGS_Model; } - /*! - * \brief Get the kind of adaptation technique. - * \return Kind of adaptation technique. - */ - unsigned short GetKind_Adaptation(void) const { return Kind_Adaptation; } - - /*! - * \brief Get the number of new elements added in the adaptation process. - * \return percentage of new elements that are going to be added in the adaptation. - */ - su2double GetNew_Elem_Adapt(void) const { return New_Elem_Adapt; } - /*! * \brief Get the kind of time integration method. * \note This is the information that the code will use, the method will @@ -5794,29 +5763,11 @@ class CConfig { */ bool GetAxisymmetric(void) const { return Axisymmetric; } - /*! - * \brief Get information about the axisymmetric frame. - * \return TRUE if there is a rotational frame; otherwise FALSE. - */ - bool GetDebugMode(void); - - /*! - * \brief Get information about there is a smoothing of the grid coordinates. - * \return TRUE if there is smoothing of the grid coordinates; otherwise FALSE. - */ - bool GetAdaptBoundary(void) const { return AdaptBoundary; } - /*! * \brief Get information about there is a smoothing of the grid coordinates. * \return TRUE if there is smoothing of the grid coordinates; otherwise FALSE. */ - bool GetSmoothNumGrid(void) const { return SmoothNumGrid; } - - /*! - * \brief Set information about there is a smoothing of the grid coordinates. - * \param[in] val_smoothnumgrid - TRUE if there is smoothing of the grid coordinates; otherwise FALSE. - */ - void SetSmoothNumGrid(bool val_smoothnumgrid) { SmoothNumGrid = val_smoothnumgrid; } + unsigned short GetSmoothNumGrid(void) const { return SmoothNumGrid; } /*! * \brief Subtract one to the index of the finest grid (full multigrid strategy). diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 647c5f58b3fc..7c9b5df7a6a2 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -67,7 +67,6 @@ enum SU2_COMPONENT { SU2_CFD = 1, /*!< \brief Running the SU2_CFD software. */ SU2_DEF = 2, /*!< \brief Running the SU2_DEF software. */ SU2_DOT = 3, /*!< \brief Running the SU2_DOT software. */ - SU2_MSH = 4, /*!< \brief Running the SU2_MSH software. */ SU2_GEO = 5, /*!< \brief Running the SU2_GEO software. */ SU2_SOL = 6 /*!< \brief Running the SU2_SOL software. */ }; @@ -1619,40 +1618,6 @@ static const MapType Sens_Map = { MakePair("SENS_AOS", SENS_AOS) }; -/*! - * \brief Types of grid adaptation/refinement - */ -enum ENUM_ADAPT { - NO_ADAPT = 0, /*!< \brief No grid adaptation. */ - FULL = 1, /*!< \brief Do a complete grid refinement of all the computational grids. */ - FULL_FLOW = 2, /*!< \brief Do a complete grid refinement of the flow grid. */ - FULL_ADJOINT = 3, /*!< \brief Do a complete grid refinement of the adjoint grid. */ - GRAD_FLOW = 5, /*!< \brief Do a gradient based grid adaptation of the flow grid. */ - GRAD_ADJOINT = 6, /*!< \brief Do a gradient based grid adaptation of the adjoint grid. */ - GRAD_FLOW_ADJ = 7, /*!< \brief Do a gradient based grid adaptation of the flow and adjoint grid. */ - COMPUTABLE = 9, /*!< \brief Apply a computable error grid adaptation. */ - REMAINING = 10, /*!< \brief Apply a remaining error grid adaptation. */ - WAKE = 12, /*!< \brief Do a grid refinement on the wake. */ - SMOOTHING = 14, /*!< \brief Do a grid smoothing of the geometry. */ - SUPERSONIC_SHOCK = 15, /*!< \brief Do a grid smoothing. */ - PERIODIC = 17 /*!< \brief Add the periodic halo cells. */ -}; -static const MapType Adapt_Map = { - MakePair("NONE", NO_ADAPT) - MakePair("FULL", FULL) - MakePair("FULL_FLOW", FULL_FLOW) - MakePair("FULL_ADJOINT", FULL_ADJOINT) - MakePair("GRAD_FLOW", GRAD_FLOW) - MakePair("GRAD_ADJOINT", GRAD_ADJOINT) - MakePair("GRAD_FLOW_ADJ", GRAD_FLOW_ADJ) - MakePair("COMPUTABLE", COMPUTABLE) - MakePair("REMAINING", REMAINING) - MakePair("WAKE", WAKE) - MakePair("SMOOTHING", SMOOTHING) - MakePair("SUPERSONIC_SHOCK", SUPERSONIC_SHOCK) - MakePair("PERIODIC", PERIODIC) -}; - /*! * \brief Types of input file formats */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 405b6f6e26b6..5d28ae540a90 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1375,7 +1375,7 @@ void CConfig::SetConfig_Options() { addStringListOption("MARKER_INTERNAL", nMarker_Internal, Marker_Internal); /* DESCRIPTION: Custom boundary marker(s) */ addStringListOption("MARKER_CUSTOM", nMarker_Custom, Marker_Custom); - /* DESCRIPTION: Periodic boundary marker(s) for use with SU2_MSH + /* DESCRIPTION: Periodic boundary marker(s) Format: ( periodic marker, donor marker, rotation_center_x, rotation_center_y, rotation_center_z, rotation_angle_x-axis, rotation_angle_y-axis, rotation_angle_z-axis, translation_x, translation_y, translation_z, ... ) */ @@ -1393,10 +1393,7 @@ void CConfig::SetConfig_Options() { /*!\brief ACTDISK_TYPE \n DESCRIPTION: Actuator Disk boundary type \n OPTIONS: see \link ActDisk_Map \endlink \n Default: VARIABLES_JUMP \ingroup Config*/ addEnumOption("ACTDISK_TYPE", Kind_ActDisk, ActDisk_Map, VARIABLES_JUMP); - /*!\brief MARKER_ACTDISK\n DESCRIPTION: Periodic boundary marker(s) for use with SU2_MSH - Format: ( periodic marker, donor marker, rotation_center_x, rotation_center_y, - rotation_center_z, rotation_angle_x-axis, rotation_angle_y-axis, - rotation_angle_z-axis, translation_x, translation_y, translation_z, ... ) \ingroup Config*/ + /*!\brief MARKER_ACTDISK\n DESCRIPTION: \ingroup Config*/ addActDiskOption("MARKER_ACTDISK", nMarker_ActDiskInlet, nMarker_ActDiskOutlet, Marker_ActDiskInlet, Marker_ActDiskOutlet, ActDisk_PressJump, ActDisk_TempJump, ActDisk_Omega); @@ -2032,21 +2029,8 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Value to move motion origins (1 or 0) */ addUShortListOption("MOVE_MOTION_ORIGIN", nMoveMotion_Origin, MoveMotion_Origin); - /*!\par CONFIG_CATEGORY: Grid adaptation \ingroup Config*/ - /*--- Options related to grid adaptation ---*/ - - /* DESCRIPTION: Kind of grid adaptation */ - addEnumOption("KIND_ADAPT", Kind_Adaptation, Adapt_Map, NO_ADAPT); - /* DESCRIPTION: Percentage of new elements (% of the original number of elements) */ - addDoubleOption("NEW_ELEMS", New_Elem_Adapt, -1.0); - /* DESCRIPTION: Scale factor for the dual volume */ - addDoubleOption("DUALVOL_POWER", DualVol_Power, 0.5); - /* DESCRIPTION: Use analytical definition for surfaces */ - addEnumOption("ANALYTICAL_SURFDEF", Analytical_Surface, Geo_Analytic_Map, NO_GEO_ANALYTIC); /* DESCRIPTION: Before each computation, implicitly smooth the nodal coordinates */ - addBoolOption("SMOOTH_GEOMETRY", SmoothNumGrid, false); - /* DESCRIPTION: Adapt the boundary elements */ - addBoolOption("ADAPT_BOUNDARY", AdaptBoundary, true); + addUnsignedShortOption("SMOOTH_GEOMETRY", SmoothNumGrid, 0); /*!\par CONFIG_CATEGORY: Aeroelastic Simulation (Typical Section Model) \ingroup Config*/ /*--- Options related to aeroelastic simulations using the Typical Section Model) ---*/ @@ -3038,7 +3022,6 @@ void CConfig::SetHeader(unsigned short val_software) const{ case SU2_CFD: cout << "| |___/\\___//___| Suite (Computational Fluid Dynamics Code) |" << endl; break; case SU2_DEF: cout << "| |___/\\___//___| Suite (Mesh Deformation Code) |" << endl; break; case SU2_DOT: cout << "| |___/\\___//___| Suite (Gradient Projection Code) |" << endl; break; - case SU2_MSH: cout << "| |___/\\___//___| Suite (Mesh Adaptation Code) |" << endl; break; case SU2_GEO: cout << "| |___/\\___//___| Suite (Geometry Definition Code) |" << endl; break; case SU2_SOL: cout << "| |___/\\___//___| Suite (Solution Exporting Code) |" << endl; break; } @@ -4758,16 +4741,6 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ Kind_ConductivityModel_Turb = NO_CONDUCTIVITY_TURB; } - /*--- Check for running SU2_MSH for periodic preprocessing, and throw - an error to report that this is no longer necessary. ---*/ - - if ((Kind_SU2 == SU2_MSH) && - (Kind_Adaptation == PERIODIC)) { - SU2_MPI::Error(string("For SU2 v7.0.0 and later, preprocessing of periodic grids by SU2_MSH\n") + - string("is no longer necessary. Please use the original mesh file (prior to SU2_MSH)\n") + - string("with the same MARKER_PERIODIC definition in the configuration file.") , CURRENT_FUNCTION); - } - /* Set a default for the size of the RECTANGLE / BOX grid sizes. */ if (nMesh_Box_Size == 0) { @@ -4980,8 +4953,7 @@ void CConfig::SetMarkers(unsigned short val_software) { int size = SINGLE_NODE; #ifdef HAVE_MPI - if (val_software != SU2_MSH) - SU2_MPI::Comm_size(MPI_COMM_WORLD, &size); + SU2_MPI::Comm_size(MPI_COMM_WORLD, &size); #endif /*--- Compute the total number of markers in the config file ---*/ @@ -5785,23 +5757,6 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) { } } - if (val_software == SU2_MSH) { - switch (Kind_Adaptation) { - case FULL: case WAKE: case FULL_FLOW: case FULL_ADJOINT: case SMOOTHING: case SUPERSONIC_SHOCK: - break; - case GRAD_FLOW: - cout << "Read flow solution from: " << Solution_FileName << "." << endl; - break; - case GRAD_ADJOINT: - cout << "Read adjoint flow solution from: " << Solution_AdjFileName << "." << endl; - break; - case GRAD_FLOW_ADJ: case COMPUTABLE: case REMAINING: - cout << "Read flow solution from: " << Solution_FileName << "." << endl; - cout << "Read adjoint flow solution from: " << Solution_AdjFileName << "." << endl; - break; - } - } - if (val_software == SU2_DEF) { cout << endl <<"---------------- Grid deformation parameters ( Zone " << iZone << " ) ----------------" << endl; cout << "Grid deformation using a linear elasticity method." << endl; @@ -6580,37 +6535,6 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) { } } - if (val_software == SU2_MSH) { - cout << endl <<"----------------- Grid adaptation strategy ( Zone " << iZone << " ) -------------------" << endl; - - switch (Kind_Adaptation) { - case NONE: break; - case PERIODIC: cout << "Grid modification to run periodic bc problems." << endl; break; - case FULL: cout << "Grid adaptation using a complete refinement." << endl; break; - case WAKE: cout << "Grid adaptation of the wake." << endl; break; - case FULL_FLOW: cout << "Flow grid adaptation using a complete refinement." << endl; break; - case FULL_ADJOINT: cout << "Adjoint grid adaptation using a complete refinement." << endl; break; - case GRAD_FLOW: cout << "Grid adaptation using gradient based strategy (density)." << endl; break; - case GRAD_ADJOINT: cout << "Grid adaptation using gradient based strategy (adjoint density)." << endl; break; - case GRAD_FLOW_ADJ: cout << "Grid adaptation using gradient based strategy (density and adjoint density)." << endl; break; - case COMPUTABLE: cout << "Grid adaptation using computable correction."<< endl; break; - case REMAINING: cout << "Grid adaptation using remaining error."<< endl; break; - case SMOOTHING: cout << "Grid smoothing using an implicit method."<< endl; break; - case SUPERSONIC_SHOCK: cout << "Grid adaptation for a supersonic shock at Mach: " << Mach <<"."<< endl; break; - } - - switch (Kind_Adaptation) { - case GRAD_FLOW: case GRAD_ADJOINT: case GRAD_FLOW_ADJ: case COMPUTABLE: case REMAINING: - cout << "Power of the dual volume in the adaptation sensor: " << DualVol_Power << endl; - cout << "Percentage of new elements in the adaptation process: " << New_Elem_Adapt << "."<< endl; - break; - } - - if (Analytical_Surface != NONE) - cout << "Use analytical definition for including points in the surfaces." << endl; - - } - cout << endl <<"-------------------- Output Information ( Zone " << iZone << " ) ----------------------" << endl; if (val_software == SU2_CFD) { @@ -6679,10 +6603,6 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) { } } - if (val_software == SU2_MSH) { - cout << "Output mesh file name: " << Mesh_Out_FileName << ". " << endl; - } - if (val_software == SU2_DOT) { if (DiscreteAdjoint) { cout << "Output Volume Sensitivity file name: " << VolSens_FileName << ". " << endl; @@ -6691,18 +6611,6 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) { cout << "Output gradient file name: " << ObjFunc_Grad_FileName << ". " << endl; } - if (val_software == SU2_MSH) { - cout << "Output mesh file name: " << Mesh_Out_FileName << ". " << endl; - cout << "Restart flow file name: " << Restart_FileName << "." << endl; - if ((Kind_Adaptation == FULL_ADJOINT) || (Kind_Adaptation == GRAD_ADJOINT) || (Kind_Adaptation == GRAD_FLOW_ADJ) || - (Kind_Adaptation == COMPUTABLE) || (Kind_Adaptation == REMAINING)) { - if (Kind_ObjFunc[0] == DRAG_COEFFICIENT) cout << "Restart adjoint file name: " << Restart_AdjFileName << "." << endl; - if (Kind_ObjFunc[0] == EQUIVALENT_AREA) cout << "Restart adjoint file name: " << Restart_AdjFileName << "." << endl; - if (Kind_ObjFunc[0] == NEARFIELD_PRESSURE) cout << "Restart adjoint file name: " << Restart_AdjFileName << "." << endl; - if (Kind_ObjFunc[0] == LIFT_COEFFICIENT) cout << "Restart adjoint file name: " << Restart_AdjFileName << "." << endl; - } - } - cout << endl <<"------------- Config File Boundary Information ( Zone " << iZone << " ) ---------------" << endl; PrintingToolbox::CTablePrinter BoundaryTable(&std::cout); diff --git a/SU2_CFD/src/output/CFlowOutput.cpp b/SU2_CFD/src/output/CFlowOutput.cpp index e7c82a93efbb..f246dd84231e 100644 --- a/SU2_CFD/src/output/CFlowOutput.cpp +++ b/SU2_CFD/src/output/CFlowOutput.cpp @@ -559,14 +559,14 @@ void CFlowOutput::SetAnalyzeSurface(CSolver *solver, CGeometry *geometry, CConfi SetHistoryOutputValue("SURFACE_STATIC_PRESSURE", Tot_Surface_Pressure); SetHistoryOutputValue("AVG_DENSITY", Tot_Surface_Density); SetHistoryOutputValue("AVG_ENTHALPY", Tot_Surface_Enthalpy); - SetHistoryOutputValue("AVG_NORMALVEL", Tot_Surface_Enthalpy); + SetHistoryOutputValue("AVG_NORMALVEL", Tot_Surface_NormalVelocity); SetHistoryOutputValue("SURFACE_UNIFORMITY", Tot_Surface_StreamVelocity2); SetHistoryOutputValue("SURFACE_SECONDARY", Tot_Surface_TransvVelocity2); SetHistoryOutputValue("SURFACE_MOM_DISTORTION", Tot_Momentum_Distortion); SetHistoryOutputValue("SURFACE_SECOND_OVER_UNIFORM", Tot_SecondOverUniformity); SetHistoryOutputValue("SURFACE_TOTAL_TEMPERATURE", Tot_Surface_TotalTemperature); SetHistoryOutputValue("SURFACE_TOTAL_PRESSURE", Tot_Surface_TotalPressure); - SetHistoryOutputValue("SURFACE_PRESSURE_DROP", Tot_Surface_PressureDrop); + SetHistoryOutputValue("SURFACE_PRESSURE_DROP", Tot_Surface_PressureDrop); if ((rank == MASTER_NODE) && !config->GetDiscrete_Adjoint() && output) { From b7b731a8971dd940dc870f12e4468e037c708e46 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 22 Jan 2021 19:15:20 +0000 Subject: [PATCH 5/9] gauss average instead of nodal (same, but different) --- .../elasticity/CFEALinearElasticity.cpp | 16 ++++++---------- .../elasticity/CFEANonlinearElasticity.cpp | 18 +++++++++++------- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp index 426200f4b598..199be05e382a 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEALinearElasticity.cpp @@ -239,8 +239,7 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, unsigned short iNode, nNode; unsigned short iDim, bDim; - /*--- Auxiliary vector ---*/ - su2double Strain[DIM_STRAIN_3D], Stress[DIM_STRAIN_3D], avgStress[DIM_STRAIN_3D] = {0.0}; + su2double avgStress[DIM_STRAIN_3D] = {0.0}; /*--- Set element properties and recompute the constitutive matrix, this is needed for multiple material cases and for correct differentiation ---*/ @@ -282,9 +281,7 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, } } - for (iVar = 0; iVar < bDim; iVar++) { - Strain[iVar] = 0.0; - } + su2double Strain[DIM_STRAIN_3D] = {0.0}; for (iNode = 0; iNode < nNode; iNode++) { @@ -319,11 +316,13 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, /*--- Compute the Stress Vector as D*epsilon ---*/ + su2double Stress[DIM_STRAIN_3D] = {0.0}; + for (iVar = 0; iVar < bDim; iVar++) { - Stress[iVar] = 0.0; for (jVar = 0; jVar < bDim; jVar++) { Stress[iVar] += D_Mat[iVar][jVar]*Strain[jVar]; } + avgStress[iVar] += Stress[iVar] / nGauss; } for (iNode = 0; iNode < nNode; iNode++) { @@ -348,10 +347,7 @@ su2double CFEALinearElasticity::Compute_Averaged_NodalStress(CElement *element, } - for (unsigned short iStress = 0; iStress < bDim; ++iStress) - for (iNode = 0; iNode < nNode; iNode++) - avgStress[iStress] += element->Get_NodalStress(iNode, iStress) / nNode; - + if (nDim == 3) std::swap(avgStress[2], avgStress[3]); auto elStress = VonMisesStress(nDim, avgStress); /*--- We only differentiate w.r.t. an avg VM stress for the element as diff --git a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp index 2dd97a983287..45f65f568fe8 100644 --- a/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp +++ b/SU2_CFD/src/numerics/elasticity/CFEANonlinearElasticity.cpp @@ -747,6 +747,8 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen unsigned short iGauss, nGauss; unsigned short iDim, iNode, nNode; + su2double avgStress[DIM_STRAIN_3D] = {0.0}; + /*--- TODO: Initialize values for the material model considered ---*/ SetElement_Properties(element, config); if (maxwell_stress) SetElectric_Properties(element, config); @@ -848,6 +850,15 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen Compute_Stress_Tensor(element, config); if (maxwell_stress) Add_MaxwellStress(element, config); + avgStress[0] += Stress_Tensor[0][0] / nGauss; + avgStress[1] += Stress_Tensor[1][1] / nGauss; + avgStress[2] += Stress_Tensor[0][1] / nGauss; + if (nDim == 3) { + avgStress[3] += Stress_Tensor[2][2] / nGauss; + avgStress[4] += Stress_Tensor[0][2] / nGauss; + avgStress[5] += Stress_Tensor[1][2] / nGauss; + } + for (iNode = 0; iNode < nNode; iNode++) { /*--- Compute the nodal stress term for each gaussian point and for each node, ---*/ @@ -879,13 +890,6 @@ su2double CFEANonlinearElasticity::Compute_Averaged_NodalStress(CElement *elemen } - su2double avgStress[DIM_STRAIN_3D] = {0.0}; - const auto nStress = (nDim == 2) ? DIM_STRAIN_2D : DIM_STRAIN_3D; - - for (unsigned short iStress = 0; iStress < nStress; ++iStress) - for (iNode = 0; iNode < nNode; iNode++) - avgStress[iStress] += element->Get_NodalStress(iNode, iStress) / nNode; - auto elStress = VonMisesStress(nDim, avgStress); /*--- We only differentiate w.r.t. an avg VM stress for the element as From 9d077df69bef2cf468d285534603ade403d99856 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 22 Jan 2021 20:28:45 +0000 Subject: [PATCH 6/9] fix p1 case --- TestCases/radiation/p1model/configp1.cfg | 1 + 1 file changed, 1 insertion(+) diff --git a/TestCases/radiation/p1model/configp1.cfg b/TestCases/radiation/p1model/configp1.cfg index 1388199e00d2..b5a60a0ffccd 100644 --- a/TestCases/radiation/p1model/configp1.cfg +++ b/TestCases/radiation/p1model/configp1.cfg @@ -23,6 +23,7 @@ INC_DENSITY_MODEL= VARIABLE INC_ENERGY_EQUATION = YES INC_DENSITY_INIT= 0.00597782417156 INC_TEMPERATURE_INIT= 288.15 +INC_VELOCITY_INIT= (0, 0, 0) INC_NONDIM = DIMENSIONAL FLUID_MODEL= INC_IDEAL_GAS From d39ec8efd61f9a36d3b4f63723af63a3fbccf8e9 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Fri, 22 Jan 2021 22:08:05 +0000 Subject: [PATCH 7/9] remove deprecated options from testcases --- .../cont_adj_rans/naca0012/turb_nasa.cfg | 19 ------------------- .../naca0012/turb_nasa_binary.cfg | 18 ------------------ .../disc_adj_euler/arina2k/Arina2KRS.cfg | 17 ----------------- .../transonic_stator_2D/transonic_stator.cfg | 6 ------ .../poiseuille/lam_poiseuille.cfg | 3 --- .../poiseuille/profile_poiseuille.cfg | 3 --- .../centrifugal_blade/centrifugal_blade.cfg | 9 --------- .../centrifugal_stage/centrifugal_stage.cfg | 10 ---------- 8 files changed, 85 deletions(-) diff --git a/TestCases/cont_adj_rans/naca0012/turb_nasa.cfg b/TestCases/cont_adj_rans/naca0012/turb_nasa.cfg index d0a39d3dc852..d9964fc5fb2b 100644 --- a/TestCases/cont_adj_rans/naca0012/turb_nasa.cfg +++ b/TestCases/cont_adj_rans/naca0012/turb_nasa.cfg @@ -257,25 +257,6 @@ CONV_CAUCHY_ELEMS= 100 % % Epsilon to control the series convergence CONV_CAUCHY_EPS= 1E-6 -% - -% ------------------------- GRID ADAPTATION STRATEGY --------------------------% -% -% Percentage of new elements (% of the original number of elements) -NEW_ELEMS= 5 -% -% Kind of grid adaptation (NONE, FULL, FULL_FLOW, GRAD_FLOW, FULL_ADJOINT, -% GRAD_ADJOINT, GRAD_FLOW_ADJ, ROBUST, -% FULL_LINEAR, COMPUTABLE, COMPUTABLE_ROBUST, -% REMAINING, WAKE, SMOOTHING, SUPERSONIC_SHOCK, -% TWOPHASE) -KIND_ADAPT= FULL_FLOW -% -% Scale factor for the dual volume -DUALVOL_POWER= 0.5 -% -% Before each computation do an implicit smoothing of the nodes coord (NO, YES) -SMOOTH_GEOMETRY= NO % ------------------------- INPUT/OUTPUT INFORMATION --------------------------% % diff --git a/TestCases/cont_adj_rans/naca0012/turb_nasa_binary.cfg b/TestCases/cont_adj_rans/naca0012/turb_nasa_binary.cfg index 662e2cda6248..66be5e33e56a 100644 --- a/TestCases/cont_adj_rans/naca0012/turb_nasa_binary.cfg +++ b/TestCases/cont_adj_rans/naca0012/turb_nasa_binary.cfg @@ -259,24 +259,6 @@ CONV_CAUCHY_ELEMS= 100 CONV_CAUCHY_EPS= 1E-6 % -% ------------------------- GRID ADAPTATION STRATEGY --------------------------% -% -% Percentage of new elements (% of the original number of elements) -NEW_ELEMS= 5 -% -% Kind of grid adaptation (NONE, FULL, FULL_FLOW, GRAD_FLOW, FULL_ADJOINT, -% GRAD_ADJOINT, GRAD_FLOW_ADJ, ROBUST, -% FULL_LINEAR, COMPUTABLE, COMPUTABLE_ROBUST, -% REMAINING, WAKE, SMOOTHING, SUPERSONIC_SHOCK, -% TWOPHASE) -KIND_ADAPT= FULL_FLOW -% -% Scale factor for the dual volume -DUALVOL_POWER= 0.5 -% -% Before each computation do an implicit smoothing of the nodes coord (NO, YES) -SMOOTH_GEOMETRY= NO - % ------------------------- INPUT/OUTPUT INFORMATION --------------------------% % % Mesh input file diff --git a/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg b/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg index 085f32cd7872..501f16e9429b 100644 --- a/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg +++ b/TestCases/disc_adj_euler/arina2k/Arina2KRS.cfg @@ -376,23 +376,6 @@ GEO_PLOT_STATIONS= NO % Geometrical evaluation mode (FUNCTION, GRADIENT) GEO_MODE= FUNCTION -% ------------------------- GRID ADAPTATION STRATEGY --------------------------% -% -% Kind of grid adaptation (NONE, PERIODIC, FULL, FULL_FLOW, GRAD_FLOW, -% FULL_ADJOINT, GRAD_ADJOINT, GRAD_FLOW_ADJ, ROBUST, -% FULL_LINEAR, COMPUTABLE, COMPUTABLE_ROBUST, -% REMAINING, WAKE, SMOOTHING, SUPERSONIC_SHOCK) -KIND_ADAPT= FULL_FLOW -% -% Percentage of new elements (% of the original number of elements) -NEW_ELEMS= 5 -% -% Scale factor for the dual volume -DUALVOL_POWER= 0.5 -% -% Adapt the boundary elements (NO, YES) -ADAPT_BOUNDARY= YES - % ----------------------- DESIGN VARIABLE PARAMETERS --------------------------% % % Kind of deformation (NO_DEFORMATION, SCALE_GRID, TRANSLATE_GRID, ROTATE_GRID, diff --git a/TestCases/disc_adj_turbomachinery/transonic_stator_2D/transonic_stator.cfg b/TestCases/disc_adj_turbomachinery/transonic_stator_2D/transonic_stator.cfg index 6e9e98762bdc..e354059f006b 100644 --- a/TestCases/disc_adj_turbomachinery/transonic_stator_2D/transonic_stator.cfg +++ b/TestCases/disc_adj_turbomachinery/transonic_stator_2D/transonic_stator.cfg @@ -187,12 +187,6 @@ MARKER_PLOTTING= (airfoil) MARKER_MONITORING= (airfoil) % % -% ------------------------- GRID ADAPTATION STRATEGY --------------------------% -% -% Kind of grid adaptation (NONE, PERIODIC) -KIND_ADAPT= PERIODIC -% -% % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % % Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES) diff --git a/TestCases/navierstokes/poiseuille/lam_poiseuille.cfg b/TestCases/navierstokes/poiseuille/lam_poiseuille.cfg index 554104dd03ba..01d165056d02 100644 --- a/TestCases/navierstokes/poiseuille/lam_poiseuille.cfg +++ b/TestCases/navierstokes/poiseuille/lam_poiseuille.cfg @@ -130,9 +130,6 @@ MARKER_PLOTTING= ( lower ) % % Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated MARKER_MONITORING= ( left, right ) -% -% Kind of adaptation (needed to create the initial periodic mesh) -KIND_ADAPT= PERIODIC % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % diff --git a/TestCases/navierstokes/poiseuille/profile_poiseuille.cfg b/TestCases/navierstokes/poiseuille/profile_poiseuille.cfg index e1704782600f..d91ca807083d 100644 --- a/TestCases/navierstokes/poiseuille/profile_poiseuille.cfg +++ b/TestCases/navierstokes/poiseuille/profile_poiseuille.cfg @@ -141,9 +141,6 @@ MARKER_PLOTTING= ( upper ) % % Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated MARKER_MONITORING= ( left, right ) -% -% Kind of adaptation (needed to create the initial periodic mesh) -KIND_ADAPT= PERIODIC % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% % diff --git a/TestCases/turbomachinery/centrifugal_blade/centrifugal_blade.cfg b/TestCases/turbomachinery/centrifugal_blade/centrifugal_blade.cfg index e33415dc7eca..caa0b59ec86d 100755 --- a/TestCases/turbomachinery/centrifugal_blade/centrifugal_blade.cfg +++ b/TestCases/turbomachinery/centrifugal_blade/centrifugal_blade.cfg @@ -157,15 +157,6 @@ MARKER_PLOTTING= ( wall1, wall2 ) MARKER_TURBO_PERFORMANCE= (inflow, outflow, BLADE) % % -% -% ------------------------- GRID ADAPTATION STRATEGY --------------------------% -% -% Kind of grid adaptation (NONE, PERIODIC) -KIND_ADAPT= PERIODIC -% -% -% -% % ----------------------- DYNAMIC MESH DEFINITION -----------------------------% % % Dynamic mesh simulation (NO, YES) diff --git a/TestCases/turbomachinery/centrifugal_stage/centrifugal_stage.cfg b/TestCases/turbomachinery/centrifugal_stage/centrifugal_stage.cfg index 40d7b49f8a36..bb572ed5c13c 100755 --- a/TestCases/turbomachinery/centrifugal_stage/centrifugal_stage.cfg +++ b/TestCases/turbomachinery/centrifugal_stage/centrifugal_stage.cfg @@ -163,16 +163,6 @@ MARKER_MONITORING= ( wall1, wall2 ) MARKER_TURBO_PERFORMANCE= (inflow, outflow, STAGE, inflow, outmix, BLADE, inmix, outflow, BLADE) % % -% -% -% ------------------------- GRID ADAPTATION STRATEGY --------------------------% -% -% Kind of grid adaptation (NONE, PERIODIC) -KIND_ADAPT= PERIODIC -% -% -% -% % ----------------------- DYNAMIC MESH DEFINITION -----------------------------% % % Dynamic mesh simulation (NO, YES) From 813ff067f109a874e8336f5e852aaab0f2cdff0b Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Mon, 25 Jan 2021 16:52:25 +0000 Subject: [PATCH 8/9] option to consider only the surface of a reference geometry --- Common/include/CConfig.hpp | 12 +++++--- Common/src/CConfig.cpp | 6 ++-- SU2_CFD/include/variables/CFEAVariable.hpp | 9 +----- SU2_CFD/include/variables/CVariable.hpp | 7 +---- SU2_CFD/src/solvers/CFEASolver.cpp | 32 +++++++++++++--------- TestCases/disc_adj_fea/configAD_fem.cfg | 4 ++- TestCases/disc_adj_fsi/configFEA.cfg | 2 ++ 7 files changed, 38 insertions(+), 34 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 3985d25e88ab..3aa76c0faa40 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -849,7 +849,7 @@ class CConfig { su2double Knowles_B, /*!< \brief Knowles material model constant B. */ Knowles_N; /*!< \brief Knowles material model constant N. */ bool DE_Effects; /*!< Application of DE effects to FE analysis */ - bool RefGeom; /*!< Read a reference geometry for optimization purposes. */ + bool RefGeom, RefGeomSurf; /*!< Read a reference geometry for optimization purposes. */ unsigned long refNodeID; /*!< \brief Global ID for the reference node (optimization). */ string RefGeom_FEMFileName; /*!< \brief File name for reference geometry. */ unsigned short RefGeom_FileFormat; /*!< \brief Mesh input format. */ @@ -2045,11 +2045,15 @@ class CConfig { su2double GetRefNode_Penalty(void) const { return RefNode_Penalty; } /*! - * \brief Decide whether it's necessary to read a reference geometry. - * \return TRUE if it's necessary to read a reference geometry, FALSE otherwise. - */ + * \brief Decide whether it's necessary to read a reference geometry. + */ bool GetRefGeom(void) const { return RefGeom; } + /*! + * \brief Consider only the surface of the reference geometry. + */ + bool GetRefGeomSurf(void) const { return RefGeomSurf; } + /*! * \brief Get the name of the file with the reference geometry of the structural problem. * \return Name of the file with the reference geometry of the structural problem. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 5d28ae540a90..45bc96adf6d6 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -2237,10 +2237,12 @@ void CConfig::SetConfig_Options() { addBoolOption("REFERENCE_GEOMETRY", RefGeom, false); /*!\brief REFERENCE_GEOMETRY_PENALTY\n DESCRIPTION: Penalty weight value for the objective function \ingroup Config*/ addDoubleOption("REFERENCE_GEOMETRY_PENALTY", RefGeom_Penalty, 1E6); - /*!\brief SOLUTION_FLOW_FILENAME \n DESCRIPTION: Restart structure input file (the file output under the filename set by RESTART_FLOW_FILENAME) \n Default: solution_flow.dat \ingroup Config */ + /*!\brief REFERENCE_GEOMETRY_FILENAME \n DESCRIPTION: Reference geometry filename \n Default: reference_geometry.dat \ingroup Config */ addStringOption("REFERENCE_GEOMETRY_FILENAME", RefGeom_FEMFileName, string("reference_geometry.dat")); - /*!\brief MESH_FORMAT \n DESCRIPTION: Mesh input file format \n OPTIONS: see \link Input_Map \endlink \n DEFAULT: SU2 \ingroup Config*/ + /*!\brief REFERENCE_GEOMETRY_FORMAT \n DESCRIPTION: Reference geometry format \n DEFAULT: SU2 \ingroup Config*/ addEnumOption("REFERENCE_GEOMETRY_FORMAT", RefGeom_FileFormat, Input_Ref_Map, SU2_REF); + /*!\brief REFERENCE_GEOMETRY_SURFACE\n DESCRIPTION: If true consider only the surfaces where loads are applied. \ingroup Config*/ + addBoolOption("REFERENCE_GEOMETRY_SURFACE", RefGeomSurf, false); /*!\brief TOTAL_DV_PENALTY\n DESCRIPTION: Penalty weight value to maintain the total sum of DV constant \ingroup Config*/ addDoubleOption("TOTAL_DV_PENALTY", DV_Penalty, 0); diff --git a/SU2_CFD/include/variables/CFEAVariable.hpp b/SU2_CFD/include/variables/CFEAVariable.hpp index c5890fad7b37..a400e09ff3a7 100644 --- a/SU2_CFD/include/variables/CFEAVariable.hpp +++ b/SU2_CFD/include/variables/CFEAVariable.hpp @@ -367,14 +367,7 @@ class CFEAVariable : public CVariable { /*! * \brief Get the pointer to the reference geometry */ - inline su2double *GetReference_Geometry(unsigned long iPoint) final { return Reference_Geometry[iPoint]; } - - /*! - * \brief Get the value of the reference geometry for the coordinate iVar - */ - inline su2double GetReference_Geometry(unsigned long iPoint, unsigned long iVar) const final { - return Reference_Geometry(iPoint,iVar); - } + inline const su2double* GetReference_Geometry(unsigned long iPoint) const final { return Reference_Geometry[iPoint]; } /*! * \brief Register the variables in the solution time_n array as input/output variable. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index 4accd4ee4580..33295cccf535 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -2308,7 +2308,7 @@ class CVariable { /*! * \brief A virtual member. */ - inline virtual su2double *GetReference_Geometry(unsigned long iPoint) {return nullptr; } + inline virtual const su2double* GetReference_Geometry(unsigned long iPoint) const { return nullptr; } /*! * \brief A virtual member. @@ -2325,11 +2325,6 @@ class CVariable { */ inline virtual su2double GetPrestretch(unsigned long iPoint, unsigned long iVar) const { return 0.0; } - /*! - * \brief A virtual member. - */ - inline virtual su2double GetReference_Geometry(unsigned long iPoint, unsigned long iVar) const { return 0.0; } - /*! * \brief A virtual member. Get the value of the undeformed coordinates. * \param[in] iDim - Index of Mesh_Coord[nDim] diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index c25fcd582ebd..8e22242e3999 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -2986,19 +2986,25 @@ void CFEASolver::Compute_OFRefGeom(CGeometry *geometry, const CConfig *config){ { su2double obj_fun_local = 0.0; - SU2_OMP_FOR_STAT(omp_chunk_size) - for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { - - for (unsigned short iVar = 0; iVar < nVar; iVar++) { - - /*--- Retrieve the value of the reference geometry ---*/ - su2double reference_geometry = nodes->GetReference_Geometry(iPoint,iVar); - - /*--- Retrieve the value of the current solution ---*/ - su2double current_solution = nodes->GetSolution(iPoint,iVar); - - /*--- The objective function is the sum of the difference between solution and difference, squared ---*/ - obj_fun_local += pow(current_solution - reference_geometry, 2); + if (!config->GetRefGeomSurf()) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + obj_fun_local += SquaredDistance(nVar, nodes->GetReference_Geometry(iPoint), nodes->GetSolution(iPoint)); + } + } + else { + for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if ((config->GetMarker_All_KindBC(iMarker) == LOAD_BOUNDARY) || + (config->GetMarker_All_KindBC(iMarker) == LOAD_DIR_BOUNDARY) || + (config->GetMarker_All_KindBC(iMarker) == FLOWLOAD_BOUNDARY)) { + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (unsigned long iVertex = 0; iVertex < geometry->GetnVertex(iMarker); ++iVertex) { + auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + + if (geometry->nodes->GetDomain(iPoint)) + obj_fun_local += SquaredDistance(nVar, nodes->GetReference_Geometry(iPoint), nodes->GetSolution(iPoint)); + } + } } } atomicAdd(obj_fun_local, objective_function); diff --git a/TestCases/disc_adj_fea/configAD_fem.cfg b/TestCases/disc_adj_fea/configAD_fem.cfg index 2877dfcbf041..b6e96e135901 100644 --- a/TestCases/disc_adj_fea/configAD_fem.cfg +++ b/TestCases/disc_adj_fea/configAD_fem.cfg @@ -4,7 +4,7 @@ % Author: R.Sanchez % % Institution: Imperial College London % % Date: 2017.11.29 % -% File Version 7.1.0 "Blackbird" % +% File Version 7.1.0 "Blackbird" % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLVER= ELASTICITY @@ -22,6 +22,8 @@ REFERENCE_GEOMETRY = YES REFERENCE_GEOMETRY_FILENAME = reference_geometry.dat REFERENCE_GEOMETRY_FORMAT = SU2 REFERENCE_GEOMETRY_PENALTY = 1E6 +% Consider only the surface +REFERENCE_GEOMETRY_SURFACE = NO READ_BINARY_RESTART=NO diff --git a/TestCases/disc_adj_fsi/configFEA.cfg b/TestCases/disc_adj_fsi/configFEA.cfg index 15d744400f95..a06a39fac287 100644 --- a/TestCases/disc_adj_fsi/configFEA.cfg +++ b/TestCases/disc_adj_fsi/configFEA.cfg @@ -24,6 +24,8 @@ ELECTRIC_FIELD_MOD = 20E5 REFERENCE_GEOMETRY = YES REFERENCE_GEOMETRY_FILENAME = reference_geometry.dat REFERENCE_GEOMETRY_FORMAT = SU2 +% Consider only the surface +REFERENCE_GEOMETRY_SURFACE = NO READ_BINARY_RESTART=NO From 30b255448b0664ac6dc9ae2e33766ee48a7c33a6 Mon Sep 17 00:00:00 2001 From: Pedro Gomes Date: Tue, 26 Jan 2021 22:20:42 +0000 Subject: [PATCH 9/9] simplify a few things --- SU2_CFD/src/solvers/CFEASolver.cpp | 204 ++++++----------------------- 1 file changed, 40 insertions(+), 164 deletions(-) diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index 8e22242e3999..84893bf9f96b 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -337,23 +337,16 @@ void CFEASolver::HybridParallelInitialization(CGeometry* geometry) { void CFEASolver::Set_ElementProperties(CGeometry *geometry, CConfig *config) { - unsigned long iElem; - unsigned long index; - unsigned long elProperties[4]; + const auto iZone = config->GetiZone(); + const auto nZone = geometry->GetnZone(); - unsigned short iZone = config->GetiZone(); - unsigned short nZone = geometry->GetnZone(); - - bool topology_mode = config->GetTopology_Optimization(); - - string filename; - ifstream properties_file; + const bool topology_mode = config->GetTopology_Optimization(); element_properties = new CProperty*[nElement]; /*--- Restart the solution from file information ---*/ - filename = config->GetFEA_FileName(); + auto filename = config->GetFEA_FileName(); /*--- If multizone, append zone name ---*/ if (nZone > 1) @@ -361,7 +354,8 @@ void CFEASolver::Set_ElementProperties(CGeometry *geometry, CConfig *config) { if (rank == MASTER_NODE) cout << "Filename: " << filename << "." << endl; - properties_file.open(filename.data(), ios::in); + ifstream properties_file; + properties_file.open(filename); /*--- In case there is no file, all elements get the same property (0) ---*/ @@ -374,7 +368,7 @@ void CFEASolver::Set_ElementProperties(CGeometry *geometry, CConfig *config) { SU2_MPI::Error("Topology mode requires an element-based properties file.",CURRENT_FUNCTION); } - for (iElem = 0; iElem < nElement; iElem++){ + for (auto iElem = 0ul; iElem < nElement; iElem++){ element_properties[iElem] = new CElementProperty(FEA_TERM, 0, 0, 0); } @@ -385,24 +379,15 @@ void CFEASolver::Set_ElementProperties(CGeometry *geometry, CConfig *config) { element_based = true; - /*--- In case this is a parallel simulation, we need to perform the - Global2Local index transformation first. ---*/ + /*--- In case this is a parallel simulation, we need to perform the Global2Local index transformation first. ---*/ - long *Global2Local = new long[geometry->GetGlobal_nElemDomain()]; + unordered_map Global2Local; - /*--- First, set all indices to a negative value by default ---*/ - - for (iElem = 0; iElem < geometry->GetGlobal_nElemDomain(); iElem++) - Global2Local[iElem] = -1; - - /*--- Now fill array with the transform values only for the points in the rank (including halos) ---*/ - - for (iElem = 0; iElem < nElement; iElem++) + for (auto iElem = 0ul; iElem < nElement; iElem++) Global2Local[geometry->elem[iElem]->GetGlobalIndex()] = iElem; /*--- Read all lines in the restart file ---*/ - long iElem_Local; unsigned long iElem_Global_Local = 0, iElem_Global = 0; string text_line; /*--- The first line is the header ---*/ @@ -420,9 +405,13 @@ void CFEASolver::Set_ElementProperties(CGeometry *geometry, CConfig *config) { Otherwise, the local index for this node on the current processor will be returned and used to instantiate the vars. ---*/ - iElem_Local = Global2Local[iElem_Global]; + auto it = Global2Local.find(iElem_Global); + + if (it != Global2Local.end()) { - if (iElem_Local >= 0) { + auto iElem_Local = it->second; + + unsigned long elProperties[4], index; if (config->GetAdvanced_FEAElementBased() || topology_mode){ point_line >> index >> elProperties[0] >> elProperties[1] >> elProperties[2] >> elProperties[3]; @@ -454,34 +443,18 @@ void CFEASolver::Set_ElementProperties(CGeometry *geometry, CConfig *config) { string("It could be empty lines at the end of the file."), CURRENT_FUNCTION); } - /*--- Close the restart file ---*/ - - properties_file.close(); - - /*--- Free memory needed for the transformation ---*/ - - delete [] Global2Local; - } } void CFEASolver::Set_Prestretch(CGeometry *geometry, CConfig *config) { - unsigned long iPoint; - unsigned long index; - - unsigned short iVar; - unsigned short iZone = config->GetiZone(); - unsigned short nZone = geometry->GetnZone(); - - string filename; - ifstream prestretch_file; - + const auto iZone = config->GetiZone(); + const auto nZone = geometry->GetnZone(); /*--- Restart the solution from file information ---*/ - filename = config->GetPrestretch_FEMFileName(); + auto filename = config->GetPrestretch_FEMFileName(); /*--- If multizone, append zone name ---*/ if (nZone > 1) @@ -489,7 +462,8 @@ void CFEASolver::Set_Prestretch(CGeometry *geometry, CConfig *config) { if (rank == MASTER_NODE) cout << "Filename: " << filename << "." << endl; - prestretch_file.open(filename.data(), ios::in); + ifstream prestretch_file; + prestretch_file.open(filename); /*--- In case there is no file ---*/ @@ -497,20 +471,15 @@ void CFEASolver::Set_Prestretch(CGeometry *geometry, CConfig *config) { SU2_MPI::Error(string("There is no FEM prestretch reference file ") + filename, CURRENT_FUNCTION); } - /*--- In case this is a parallel simulation, we need to perform the - Global2Local index transformation first. ---*/ - - map Global2Local; - map::const_iterator MI; + /*--- Make a global to local map that also covers halo nodes (the one in geometry does not). ---*/ - /*--- Now fill array with the transform values only for local points ---*/ + unordered_map Global2Local; - for (iPoint = 0; iPoint < nPointDomain; iPoint++) + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) Global2Local[geometry->nodes->GetGlobalIndex(iPoint)] = iPoint; /*--- Read all lines in the restart file ---*/ - long iPoint_Local; unsigned long iPoint_Global_Local = 0, iPoint_Global = 0; string text_line; /*--- The first line is the header ---*/ @@ -523,17 +492,20 @@ void CFEASolver::Set_Prestretch(CGeometry *geometry, CConfig *config) { /*--- Retrieve local index. If this node from the restart file lives on the current processor, we will load and instantiate the vars. ---*/ - MI = Global2Local.find(iPoint_Global); - if (MI != Global2Local.end()) { + auto it = Global2Local.find(iPoint_Global); + + if (it != Global2Local.end()) { - iPoint_Local = Global2Local[iPoint_Global]; + auto iPoint_Local = it->second; su2double Sol[MAXNVAR] = {0.0}; + unsigned long index; if (nDim == 2) point_line >> Sol[0] >> Sol[1] >> index; if (nDim == 3) point_line >> Sol[0] >> Sol[1] >> Sol[2] >> index; - for (iVar = 0; iVar < nVar; iVar++) nodes->SetPrestretch(iPoint_Local,iVar, Sol[iVar]); + for (unsigned short iVar = 0; iVar < nVar; iVar++) + nodes->SetPrestretch(iPoint_Local, iVar, Sol[iVar]); iPoint_Global_Local++; } @@ -542,99 +514,28 @@ void CFEASolver::Set_Prestretch(CGeometry *geometry, CConfig *config) { /*--- Detect a wrong solution file ---*/ - if (iPoint_Global_Local != nPointDomain) { + if (iPoint_Global_Local != nPoint) { SU2_MPI::Error(string("The solution file ") + filename + string(" doesn't match with the mesh file!\n") + string("It could be empty lines at the end of the file."), CURRENT_FUNCTION); } - /*--- Close the restart file ---*/ - - prestretch_file.close(); - -#ifdef HAVE_MPI - /*--- We need to communicate here the prestretched geometry for the halo nodes. ---*/ - /*--- We avoid creating a new function as this may be reformatted. ---*/ - - unsigned short iMarker, MarkerS, MarkerR; - unsigned long iVertex, nVertexS, nVertexR, nBufferS_Vector, nBufferR_Vector; - su2double *Buffer_Receive_U = nullptr, *Buffer_Send_U = nullptr; - - int send_to, receive_from; - - for (iMarker = 0; iMarker < nMarker; iMarker++) { - - if ((config->GetMarker_All_KindBC(iMarker) == SEND_RECEIVE) && - (config->GetMarker_All_SendRecv(iMarker) > 0)) { - - MarkerS = iMarker; MarkerR = iMarker+1; - - send_to = config->GetMarker_All_SendRecv(MarkerS)-1; - receive_from = abs(config->GetMarker_All_SendRecv(MarkerR))-1; - - nVertexS = geometry->nVertex[MarkerS]; nVertexR = geometry->nVertex[MarkerR]; - nBufferS_Vector = nVertexS*nVar; nBufferR_Vector = nVertexR*nVar; - - /*--- Allocate Receive and send buffers ---*/ - Buffer_Receive_U = new su2double [nBufferR_Vector]; - Buffer_Send_U = new su2double[nBufferS_Vector]; - - /*--- Copy the solution that should be sent ---*/ - for (iVertex = 0; iVertex < nVertexS; iVertex++) { - iPoint = geometry->vertex[MarkerS][iVertex]->GetNode(); - for (iVar = 0; iVar < nVar; iVar++) - Buffer_Send_U[iVar*nVertexS+iVertex] = nodes->GetPrestretch(iPoint,iVar); - } - - /*--- Send/Receive information using Sendrecv ---*/ - SU2_MPI::Sendrecv(Buffer_Send_U, nBufferS_Vector, MPI_DOUBLE, send_to, 0, - Buffer_Receive_U, nBufferR_Vector, MPI_DOUBLE, receive_from, 0, MPI_COMM_WORLD, nullptr); - - /*--- Deallocate send buffer ---*/ - delete [] Buffer_Send_U; - - /*--- Do the coordinate transformation ---*/ - for (iVertex = 0; iVertex < nVertexR; iVertex++) { - - /*--- Find point and its type of transformation ---*/ - iPoint = geometry->vertex[MarkerR][iVertex]->GetNode(); - - /*--- Store received values back into the variable. ---*/ - for (iVar = 0; iVar < nVar; iVar++) - nodes->SetPrestretch(iPoint, iVar, Buffer_Receive_U[iVar*nVertexR+iVertex]); - - } - - /*--- Deallocate receive buffer ---*/ - delete [] Buffer_Receive_U; - - } - - } -#endif - } void CFEASolver::Set_ReferenceGeometry(CGeometry *geometry, CConfig *config) { - unsigned long iPoint; - - unsigned short iVar; - unsigned short iZone = config->GetiZone(); - unsigned short file_format = config->GetRefGeom_FileFormat(); - - string filename; - ifstream reference_file; - + const auto iZone = config->GetiZone(); + const auto file_format = config->GetRefGeom_FileFormat(); /*--- Restart the solution from file information ---*/ - filename = config->GetRefGeom_FEMFileName(); + auto filename = config->GetRefGeom_FEMFileName(); /*--- If multizone, append zone name ---*/ filename = config->GetMultizone_FileName(filename, iZone, ".csv"); - reference_file.open(filename.data(), ios::in); + ifstream reference_file; + reference_file.open(filename); /*--- In case there is no file ---*/ @@ -644,24 +545,8 @@ void CFEASolver::Set_ReferenceGeometry(CGeometry *geometry, CConfig *config) { if (rank == MASTER_NODE) cout << "Filename: " << filename << " and format " << file_format << "." << endl; - /*--- In case this is a parallel simulation, we need to perform the - Global2Local index transformation first. ---*/ - - long *Global2Local = new long[geometry->GetGlobal_nPointDomain()]; - - /*--- First, set all indices to a negative value by default ---*/ - - for (iPoint = 0; iPoint < geometry->GetGlobal_nPointDomain(); iPoint++) - Global2Local[iPoint] = -1; - - /*--- Now fill array with the transform values only for local points ---*/ - - for (iPoint = 0; iPoint < nPointDomain; iPoint++) - Global2Local[geometry->nodes->GetGlobalIndex(iPoint)] = iPoint; - /*--- Read all lines in the restart file ---*/ - long iPoint_Local; unsigned long iPoint_Global_Local = 0, iPoint_Global = 0; string text_line; /*--- The first line is the header ---*/ @@ -677,7 +562,7 @@ void CFEASolver::Set_ReferenceGeometry(CGeometry *geometry, CConfig *config) { Otherwise, the local index for this node on the current processor will be returned and used to instantiate the vars. ---*/ - iPoint_Local = Global2Local[iPoint_Global]; + auto iPoint_Local = geometry->GetGlobal_to_Local_Point(iPoint_Global); if (iPoint_Local >= 0) { @@ -692,7 +577,7 @@ void CFEASolver::Set_ReferenceGeometry(CGeometry *geometry, CConfig *config) { Sol[2] = PrintingToolbox::stod(point_line[6]); } - for (iVar = 0; iVar < nVar; iVar++) + for (unsigned short iVar = 0; iVar < nVar; iVar++) nodes->SetReference_Geometry(iPoint_Local, iVar, Sol[iVar]); iPoint_Global_Local++; @@ -707,14 +592,6 @@ void CFEASolver::Set_ReferenceGeometry(CGeometry *geometry, CConfig *config) { string("It could be empty lines at the end of the file."), CURRENT_FUNCTION); } - /*--- Close the restart file ---*/ - - reference_file.close(); - - /*--- Free memory needed for the transformation ---*/ - - delete [] Global2Local; - } void CFEASolver::Set_VertexEliminationSchedule(CGeometry *geometry, const vector& markers) { @@ -3387,9 +3264,8 @@ void CFEASolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config) if (rank == MASTER_NODE) { string filename = config->GetTopology_Optim_FileName(); ofstream file; - file.open(filename.c_str()); + file.open(filename); for(iElem=0; iElem