From c5e97f862bb0bdd93fd48f9f666ab2993724866f Mon Sep 17 00:00:00 2001 From: Pedro Gomes <38071223+pcarruscag@users.noreply.github.com> Date: Thu, 24 Aug 2023 08:50:19 -0700 Subject: [PATCH] Multizone py wrapper example (FSI) (#2111) * checkpoint naming and update some functions * add FSI example with additional custom load via the python wrapper * remove mesh files * fix segfault * add some arm64 values * disable polar naca on arm64 --- SU2_CFD/include/drivers/CDriverBase.hpp | 21 +++-- SU2_CFD/include/drivers/CMultizoneDriver.hpp | 20 ++++- .../interfaces/fsi/CFlowTractionInterface.hpp | 5 +- .../include/variables/CFEABoundVariable.hpp | 5 -- SU2_CFD/include/variables/CVariable.hpp | 5 -- SU2_CFD/src/drivers/CDriverBase.cpp | 10 +-- SU2_CFD/src/drivers/CMultizoneDriver.cpp | 43 +++++----- SU2_CFD/src/drivers/CSinglezoneDriver.cpp | 15 ++-- .../fsi/CDiscAdjFlowTractionInterface.cpp | 6 +- .../interfaces/fsi/CFlowTractionInterface.cpp | 28 +++++-- SU2_CFD/src/output/CElasticityOutput.cpp | 6 +- SU2_CFD/src/python_wrapper_structure.cpp | 60 +++++++------- SU2_CFD/src/variables/CFEABoundVariable.cpp | 2 - TestCases/parallel_regression.py | 13 +++ TestCases/py_wrapper/dyn_fsi/config.cfg | 31 +++++++ TestCases/py_wrapper/dyn_fsi/configFEA.cfg | 54 ++++++++++++ TestCases/py_wrapper/dyn_fsi/configFlow.cfg | 80 ++++++++++++++++++ TestCases/py_wrapper/dyn_fsi/run.py | 83 +++++++++++++++++++ TestCases/serial_regression.py | 2 + TestCases/vandv.py | 1 + 20 files changed, 386 insertions(+), 104 deletions(-) create mode 100644 TestCases/py_wrapper/dyn_fsi/config.cfg create mode 100644 TestCases/py_wrapper/dyn_fsi/configFEA.cfg create mode 100644 TestCases/py_wrapper/dyn_fsi/configFlow.cfg create mode 100644 TestCases/py_wrapper/dyn_fsi/run.py diff --git a/SU2_CFD/include/drivers/CDriverBase.hpp b/SU2_CFD/include/drivers/CDriverBase.hpp index 6ee8a860ca5..01c3ae7a58d 100644 --- a/SU2_CFD/include/drivers/CDriverBase.hpp +++ b/SU2_CFD/include/drivers/CDriverBase.hpp @@ -54,7 +54,7 @@ class CDriverBase { UsedTime; /*!< \brief Elapsed time between Start and Stop point of the timer. */ unsigned long TimeIter; - unsigned short selected_iZone = ZONE_0; /*!< \brief Selected zone for the driver. Defaults to ZONE_0 */ + unsigned short selected_zone = ZONE_0; /*!< \brief Selected zone for the driver. Defaults to ZONE_0 */ unsigned short iMesh, /*!< \brief Iterator on mesh levels. */ iZone, /*!< \brief Iterator on zones. */ nZone, /*!< \brief Total number of zones in the problem. */ @@ -227,7 +227,7 @@ class CDriverBase { SU2_MPI::Error("Initial coordinates are only available with DEFORM_MESH= YES", CURRENT_FUNCTION); } auto* coords = - const_cast(solver_container[selected_iZone][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord()); + const_cast(solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord()); return CPyWrapperMatrixView(*coords, "InitialCoordinates", true); } @@ -241,7 +241,7 @@ class CDriverBase { if (iMarker >= GetNumberMarkers()) SU2_MPI::Error("Marker index exceeds size.", CURRENT_FUNCTION); auto* coords = - const_cast(solver_container[selected_iZone][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord()); + const_cast(solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->GetNodes()->GetMesh_Coord()); return CPyWrapperMarkerMatrixView(*coords, main_geometry->vertex[iMarker], main_geometry->GetnVertex(iMarker), "MarkerInitialCoordinates", true); } @@ -550,16 +550,21 @@ class CDriverBase { } /*! - * \brief Selects zone to be used for Driver operation + * \brief Selects zone to be used for python driver operations. * \param[in] iZone - Zone identifier. */ inline void SelectZone(unsigned short iZone) { if (iZone >= nZone) SU2_MPI::Error("Zone index out of range", CURRENT_FUNCTION); - selected_iZone = iZone; - main_geometry = geometry_container[selected_iZone][INST_0][MESH_0]; - main_config = config_container[selected_iZone]; + selected_zone = iZone; + main_geometry = geometry_container[selected_zone][INST_0][MESH_0]; + main_config = config_container[selected_zone]; } + /*! + * \brief Returns the index of the zone selected for python driver operations. + */ + inline unsigned short SelectedZone() const { return selected_zone; } + /*! * \brief Get the wall normal heat flux at a vertex on a specified marker of the flow or heat solver. * \note This can be the output of a heat or flow solver in a CHT setting. @@ -707,7 +712,7 @@ class CDriverBase { if (iMarker < std::numeric_limits::max() && iMarker > GetNumberMarkers()) { SU2_MPI::Error("Marker index exceeds size.", CURRENT_FUNCTION); } - auto* solver = solver_container[selected_iZone][INST_0][MESH_0][iSolver]; + auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSolver]; if (solver == nullptr) SU2_MPI::Error("The selected solver does not exist.", CURRENT_FUNCTION); return solver; } diff --git a/SU2_CFD/include/drivers/CMultizoneDriver.hpp b/SU2_CFD/include/drivers/CMultizoneDriver.hpp index afff04b8eac..adf9bf9d617 100644 --- a/SU2_CFD/include/drivers/CMultizoneDriver.hpp +++ b/SU2_CFD/include/drivers/CMultizoneDriver.hpp @@ -122,7 +122,7 @@ class CMultizoneDriver : public CDriver { /*! * \brief Destructor of the class. */ - ~CMultizoneDriver(void) override; + ~CMultizoneDriver() override; /*! * \brief [Overload] Launch the computation for multizone problems. @@ -130,10 +130,26 @@ class CMultizoneDriver : public CDriver { void StartSolver() override; /*! - * \brief Preprocess the multizone iteration + * \brief Preprocess the multizone iteration. */ void Preprocess(unsigned long TimeIter) override; + /*! + * \brief Solves one time iteration. + */ + void Run() override { + switch (driver_config->GetKind_MZSolver()){ + case ENUM_MULTIZONE::MZ_BLOCK_GAUSS_SEIDEL: RunGaussSeidel(); break; // Block Gauss-Seidel iteration + case ENUM_MULTIZONE::MZ_BLOCK_JACOBI: RunJacobi(); break; // Block-Jacobi iteration + } + } + + /*! + * \brief Placeholder for post processing operations to make the interface + * of this driver identical to CSinglezoneDriver. + */ + void Postprocess() {} + /*! * \brief Update the dual-time solution within multiple zones. */ diff --git a/SU2_CFD/include/interfaces/fsi/CFlowTractionInterface.hpp b/SU2_CFD/include/interfaces/fsi/CFlowTractionInterface.hpp index 2074319c110..d3ea4a6a59d 100644 --- a/SU2_CFD/include/interfaces/fsi/CFlowTractionInterface.hpp +++ b/SU2_CFD/include/interfaces/fsi/CFlowTractionInterface.hpp @@ -43,8 +43,11 @@ class CFlowTractionInterface : public CInterface { /*! * \brief Sets the dimensional factor for pressure and the consistent_interpolation flag. * \param[in] flow_config - Definition of the fluid (donor) problem. + * \param[in] struct_config - Definition of the structural (target) problem. + * \param[in] geometry - FEA geometry. + * \param[in] solution - FEA solver. */ - void Preprocess(const CConfig *flow_config); + void Preprocess(const CConfig *flow_config, const CConfig *struct_config, CGeometry *geometry, CSolver *solution); /*! * \brief Computes vertex areas (FEA side) for when tractions need to be integrated. diff --git a/SU2_CFD/include/variables/CFEABoundVariable.hpp b/SU2_CFD/include/variables/CFEABoundVariable.hpp index e0f8c273cc1..cdac1cdd476 100644 --- a/SU2_CFD/include/variables/CFEABoundVariable.hpp +++ b/SU2_CFD/include/variables/CFEABoundVariable.hpp @@ -149,11 +149,6 @@ class CFEABoundVariable final : public CFEAVariable { return FlowTraction_n(iPoint,iVar); } - /*! - * \brief Clear the flow traction residual - */ - void Clear_FlowTraction() override; - /*! * \brief Register the flow tractions as input variable. */ diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index ec2557c8088..af0d70f2673 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -1508,11 +1508,6 @@ class CVariable { */ inline virtual su2double Get_FlowTraction_n(unsigned long iPoint, unsigned long iVar) const { return 0.0; } - /*! - * \brief A virtual member. - */ - inline virtual void Clear_FlowTraction() {} - /*! * \brief A virtual member. */ diff --git a/SU2_CFD/src/drivers/CDriverBase.cpp b/SU2_CFD/src/drivers/CDriverBase.cpp index 603e2f734f8..82bb3c3a477 100644 --- a/SU2_CFD/src/drivers/CDriverBase.cpp +++ b/SU2_CFD/src/drivers/CDriverBase.cpp @@ -391,14 +391,14 @@ vector CDriverBase::GetMarkerVertexNormals(unsigned short iMarker } void CDriverBase::CommunicateMeshDisplacements() { - solver_container[selected_iZone][INST_0][MESH_0][MESH_SOL]->InitiateComms(main_geometry, main_config, MESH_DISPLACEMENTS); - solver_container[selected_iZone][INST_0][MESH_0][MESH_SOL]->CompleteComms(main_geometry, main_config, MESH_DISPLACEMENTS); + solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->InitiateComms(main_geometry, main_config, MESH_DISPLACEMENTS); + solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->CompleteComms(main_geometry, main_config, MESH_DISPLACEMENTS); } map CDriverBase::GetSolverIndices() const { map indexMap; for (auto iSol = 0u; iSol < MAX_SOLS; iSol++) { - const auto* solver = solver_container[selected_iZone][INST_0][MESH_0][iSol]; + const auto* solver = solver_container[selected_zone][INST_0][MESH_0][iSol]; if (solver != nullptr) { if (solver->GetSolverName().empty()) SU2_MPI::Error("Solver name was not defined.", CURRENT_FUNCTION); indexMap[solver->GetSolverName()] = iSol; @@ -408,7 +408,7 @@ map CDriverBase::GetSolverIndices() const { } std::map CDriverBase::GetFEASolutionIndices() const { - if (solver_container[selected_iZone][INST_0][MESH_0][FEA_SOL] == nullptr) { + if (solver_container[selected_zone][INST_0][MESH_0][FEA_SOL] == nullptr) { SU2_MPI::Error("The FEA solver does not exist.", CURRENT_FUNCTION); } const auto nDim = main_geometry->GetnDim(); @@ -429,7 +429,7 @@ std::map CDriverBase::GetFEASolutionIndices() const { } map CDriverBase::GetPrimitiveIndices() const { - if (solver_container[selected_iZone][INST_0][MESH_0][FLOW_SOL] == nullptr) { + if (solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL] == nullptr) { SU2_MPI::Error("The flow solver does not exist.", CURRENT_FUNCTION); } return PrimitiveNameToIndexMap(CPrimitiveIndices( diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index b3a7af5e369..14c31aa9ea3 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -182,10 +182,7 @@ void CMultizoneDriver::StartSolver() { /*--- Run a block iteration of the multizone problem. ---*/ - switch (driver_config->GetKind_MZSolver()){ - case ENUM_MULTIZONE::MZ_BLOCK_GAUSS_SEIDEL: RunGaussSeidel(); break; // Block Gauss-Seidel iteration - case ENUM_MULTIZONE::MZ_BLOCK_JACOBI: RunJacobi(); break; // Block-Jacobi iteration - } + Run(); /*--- Update the solution for dual time stepping strategy ---*/ @@ -646,11 +643,11 @@ void CMultizoneDriver::SetTurboPerformance() { } } -bool CMultizoneDriver::Monitor(unsigned long TimeIter){ +bool CMultizoneDriver::Monitor(unsigned long TimeIter) { /*--- Check whether the inner solver has converged --- */ - if (driver_config->GetTime_Domain() == NO){ + if (driver_config->GetTime_Domain() == NO) { const auto OuterIter = driver_config->GetOuterIter(); const auto nOuterIter = driver_config->GetnOuter_Iter(); @@ -670,29 +667,27 @@ bool CMultizoneDriver::Monitor(unsigned long TimeIter){ } // i.e. unsteady simulation - /*--- Check whether the outer time integration has reached the final time ---*/ - const auto TimeConvergence = GetTimeConvergence(); + /*--- Check whether the outer time integration has reached the final time. ---*/ + const auto TimeConvergence = GetTimeConvergence(); - const auto nTimeIter = driver_config->GetnTime_Iter(); - const auto MaxTime = driver_config->GetMax_Time(); - const auto CurTime = driver_output->GetHistoryFieldValue("CUR_TIME"); + const auto nTimeIter = driver_config->GetnTime_Iter(); + const auto MaxTime = driver_config->GetMax_Time(); + const auto CurTime = driver_output->GetHistoryFieldValue("CUR_TIME"); - const bool FinalTimeReached = (CurTime >= MaxTime); - const bool MaxIterationsReached = (TimeIter+1 >= nTimeIter); + const bool FinalTimeReached = (CurTime >= MaxTime); + const bool MaxIterationsReached = (TimeIter+1 >= nTimeIter); - if ((TimeConvergence || FinalTimeReached || MaxIterationsReached) && (rank == MASTER_NODE)){ - cout << "\n----------------------------- Solver Exit -------------------------------"; - if (TimeConvergence) cout << "\nAll windowed time-averaged convergence criteria are fullfilled." << endl; - if (FinalTimeReached) cout << "\nMaximum time reached (MAX_TIME = " << MaxTime << "s)." << endl; - else cout << "\nMaximum number of time iterations reached (TIME_ITER = " << nTimeIter << ")." << endl; - cout << "-------------------------------------------------------------------------" << endl; - } - - return (FinalTimeReached || MaxIterationsReached); - + if ((TimeConvergence || FinalTimeReached || MaxIterationsReached) && (rank == MASTER_NODE)){ + cout << "\n----------------------------- Solver Exit -------------------------------"; + if (TimeConvergence) cout << "\nAll windowed time-averaged convergence criteria are fullfilled." << endl; + if (FinalTimeReached) cout << "\nMaximum time reached (MAX_TIME = " << MaxTime << "s)." << endl; + else cout << "\nMaximum number of time iterations reached (TIME_ITER = " << nTimeIter << ")." << endl; + cout << "-------------------------------------------------------------------------" << endl; + } - if (rank == MASTER_NODE) SetTurboPerformance(); + if (rank == MASTER_NODE && driver_config->GetBoolTurbomachinery()) SetTurboPerformance(); + return (FinalTimeReached || MaxIterationsReached); } bool CMultizoneDriver::GetTimeConvergence() const{ diff --git a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp index c1f2c2519c8..f9fa11e87d6 100644 --- a/SU2_CFD/src/drivers/CSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CSinglezoneDriver.cpp @@ -92,12 +92,6 @@ void CSinglezoneDriver::StartSolver() { Output(TimeIter); - /*--- Save iteration solution for libROM ---*/ - if (config_container[MESH_0]->GetSave_libROM()) { - solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SavelibROM(geometry_container[ZONE_0][INST_0][MESH_0], - config_container[ZONE_0], StopCalc); - } - /*--- If the convergence criteria has been met, terminate the simulation. ---*/ if (StopCalc) break; @@ -208,7 +202,14 @@ void CSinglezoneDriver::Output(unsigned long TimeIter) { solver_container[ZONE_0][INST_0][MESH_0], TimeIter, StopCalc); - if (wrote_files){ + /*--- Save iteration solution for libROM ---*/ + if (config_container[MESH_0]->GetSave_libROM()) { + solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->SavelibROM(geometry_container[ZONE_0][INST_0][MESH_0], + config_container[ZONE_0], StopCalc); + wrote_files = true; + } + + if (wrote_files) { StopTime = SU2_MPI::Wtime(); diff --git a/SU2_CFD/src/interfaces/fsi/CDiscAdjFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CDiscAdjFlowTractionInterface.cpp index b8a48ae8fc4..2db45dd5825 100644 --- a/SU2_CFD/src/interfaces/fsi/CDiscAdjFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CDiscAdjFlowTractionInterface.cpp @@ -40,11 +40,7 @@ void CDiscAdjFlowTractionInterface::GetPhysical_Constants(CSolver *flow_solution CGeometry *flow_geometry, CGeometry *struct_geometry, const CConfig *flow_config, const CConfig *struct_config){ - /*--- We have to clear the traction before applying it, because we are "adding" to node and not "setting" ---*/ - - struct_solution->GetNodes()->Clear_FlowTraction(); - - Preprocess(flow_config); + Preprocess(flow_config, struct_config, struct_geometry, struct_solution); if (!conservative) ComputeVertexAreas(struct_config, struct_geometry, struct_solution); diff --git a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp index 1aaa0d422fe..88ee4a5e41f 100644 --- a/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CFlowTractionInterface.cpp @@ -31,6 +31,7 @@ #include "../../../../Common/include/geometry/CGeometry.hpp" #include "../../../include/solvers/CSolver.hpp" #include "../../../../Common/include/toolboxes/geometry_toolbox.hpp" +#include "../../../Common/include/interface_interpolation/CInterpolator.hpp" #include CFlowTractionInterface::CFlowTractionInterface(unsigned short val_nVar, unsigned short val_nConst, @@ -39,7 +40,26 @@ CFlowTractionInterface::CFlowTractionInterface(unsigned short val_nVar, unsigned conservative(conservative_) { } -void CFlowTractionInterface::Preprocess(const CConfig *flow_config) { +void CFlowTractionInterface::Preprocess(const CConfig *flow_config, const CConfig *struct_config, + CGeometry *struct_geometry, CSolver *struct_solution) { + + /*--- Clear the tractions only on the markers involved in interface, fluid tractions + * on other markers can be specified via e.g. the python wrapper. ---*/ + + for (auto iMarkerInt = 0u; iMarkerInt < struct_config->GetMarker_n_ZoneInterface()/2; iMarkerInt++) { + const auto markDonor = flow_config->FindInterfaceMarker(iMarkerInt); + const auto markTarget = struct_config->FindInterfaceMarker(iMarkerInt); + + if(!CInterpolator::CheckInterfaceBoundary(markDonor, markTarget)) continue; + if (markTarget < 0) continue; + + for (auto iVertex = 0ul; iVertex < struct_geometry->GetnVertex(markTarget); iVertex++) { + const auto iPoint = struct_geometry->vertex[markTarget][iVertex]->GetNode(); + if (!struct_geometry->nodes->GetDomain(iPoint)) continue; + su2double zeros[3] = {}; + struct_solution->GetNodes()->Set_FlowTraction(iPoint, zeros); + } + } /*--- Compute the constant factor to dimensionalize pressure and shear stress. ---*/ const su2double *Velocity_ND, *Velocity_Real; @@ -121,11 +141,7 @@ void CFlowTractionInterface::GetPhysical_Constants(CSolver *flow_solution, CSolv CGeometry *flow_geometry, CGeometry *struct_geometry, const CConfig *flow_config, const CConfig *struct_config) { - /*--- We have to clear the traction before applying it, because we are "adding" to node and not "setting" ---*/ - - struct_solution->GetNodes()->Clear_FlowTraction(); - - Preprocess(flow_config); + Preprocess(flow_config, struct_config, struct_geometry, struct_solution); if (!conservative) ComputeVertexAreas(struct_config, struct_geometry, struct_solution); diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index fe5efe5b390..71726c6ebb3 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -98,8 +98,10 @@ CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COu /*--- Set the default convergence field --- */ - if (convFields.empty() ) convFields.emplace_back("RMS_DISP_X"); - + if (convFields.empty()) { + if (linear_analysis) convFields.emplace_back("RMS_DISP_X"); + if (nonlinear_analysis) convFields.emplace_back("RMS_UTOL"); + } } void CElasticityOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolver **solver) { diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index a189217aa63..ac1b774e808 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -62,15 +62,15 @@ void CDriver::PreprocessPythonInterface(CConfig** config, CGeometry**** geometry /* Functions to obtain global parameters from SU2 (time steps, delta t, etc.) */ ////////////////////////////////////////////////////////////////////////////////// -unsigned long CDriver::GetNumberTimeIter() const { return config_container[selected_iZone]->GetnTime_Iter(); } +unsigned long CDriver::GetNumberTimeIter() const { return config_container[selected_zone]->GetnTime_Iter(); } unsigned long CDriver::GetTimeIter() const { return TimeIter; } passivedouble CDriver::GetUnsteadyTimeStep() const { - return SU2_TYPE::GetValue(config_container[selected_iZone]->GetTime_Step()); + return SU2_TYPE::GetValue(config_container[selected_zone]->GetTime_Step()); } -string CDriver::GetSurfaceFileName() const { return config_container[selected_iZone]->GetSurfCoeff_FileName(); } +string CDriver::GetSurfaceFileName() const { return config_container[selected_zone]->GetSurfCoeff_FileName(); } //////////////////////////////////////////////////////////////////////////////// /* Functions related to the management of markers */ @@ -78,12 +78,12 @@ string CDriver::GetSurfaceFileName() const { return config_container[selected_iZ void CDriver::SetHeatSourcePosition(passivedouble alpha, passivedouble pos_x, passivedouble pos_y, passivedouble pos_z) { - CSolver* solver = solver_container[selected_iZone][INST_0][MESH_0][RAD_SOL]; + CSolver* solver = solver_container[selected_zone][INST_0][MESH_0][RAD_SOL]; - config_container[selected_iZone]->SetHeatSource_Rot_Z(alpha); - config_container[selected_iZone]->SetHeatSource_Center(pos_x, pos_y, pos_z); + config_container[selected_zone]->SetHeatSource_Rot_Z(alpha); + config_container[selected_zone]->SetHeatSource_Center(pos_x, pos_y, pos_z); - solver->SetVolumetricHeatSource(geometry_container[selected_iZone][INST_0][MESH_0], config_container[selected_iZone]); + solver->SetVolumetricHeatSource(geometry_container[selected_zone][INST_0][MESH_0], config_container[selected_zone]); } void CDriver::SetInletAngle(unsigned short iMarker, passivedouble alpha) { @@ -91,20 +91,20 @@ void CDriver::SetInletAngle(unsigned short iMarker, passivedouble alpha) { unsigned long iVertex; - for (iVertex = 0; iVertex < geometry_container[selected_iZone][INST_0][MESH_0]->nVertex[iMarker]; iVertex++) { - solver_container[selected_iZone][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, 0, cos(alpha_rad)); - solver_container[selected_iZone][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, 1, sin(alpha_rad)); + for (iVertex = 0; iVertex < geometry_container[selected_zone][INST_0][MESH_0]->nVertex[iMarker]; iVertex++) { + solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, 0, cos(alpha_rad)); + solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL]->SetInlet_FlowDir(iMarker, iVertex, 1, sin(alpha_rad)); } } void CDriver::SetFarFieldAoA(const passivedouble AoA) { - config_container[ZONE_0]->SetAoA(AoA); - solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->UpdateFarfieldVelocity(config_container[ZONE_0]); + config_container[selected_zone]->SetAoA(AoA); + solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL]->UpdateFarfieldVelocity(config_container[selected_zone]); } void CDriver::SetFarFieldAoS(const passivedouble AoS) { - config_container[ZONE_0]->SetAoS(AoS); - solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->UpdateFarfieldVelocity(config_container[ZONE_0]); + config_container[selected_zone]->SetAoS(AoS); + solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL]->UpdateFarfieldVelocity(config_container[selected_zone]); } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -116,22 +116,22 @@ void CSinglezoneDriver::SetInitialMesh() { SU2_OMP_PARALLEL { for (iMesh = 0u; iMesh <= main_config->GetnMGLevels(); iMesh++) { - SU2_OMP_FOR_STAT(roundUpDiv(geometry_container[selected_iZone][INST_0][iMesh]->GetnPoint(), omp_get_max_threads())) - for (auto iPoint = 0ul; iPoint < geometry_container[selected_iZone][INST_0][iMesh]->GetnPoint(); iPoint++) { + SU2_OMP_FOR_STAT(roundUpDiv(geometry_container[selected_zone][INST_0][iMesh]->GetnPoint(), omp_get_max_threads())) + for (auto iPoint = 0ul; iPoint < geometry_container[selected_zone][INST_0][iMesh]->GetnPoint(); iPoint++) { /*--- Overwrite fictitious velocities. ---*/ su2double Grid_Vel[3] = {0.0, 0.0, 0.0}; /*--- Set the grid velocity for this coarse node. ---*/ - geometry_container[selected_iZone][INST_0][iMesh]->nodes->SetGridVel(iPoint, Grid_Vel); + geometry_container[selected_zone][INST_0][iMesh]->nodes->SetGridVel(iPoint, Grid_Vel); } END_SU2_OMP_FOR /*--- Push back the volume. ---*/ - geometry_container[selected_iZone][INST_0][iMesh]->nodes->SetVolume_n(); - geometry_container[selected_iZone][INST_0][iMesh]->nodes->SetVolume_nM1(); + geometry_container[selected_zone][INST_0][iMesh]->nodes->SetVolume_n(); + geometry_container[selected_zone][INST_0][iMesh]->nodes->SetVolume_nM1(); } /*--- Push back the solution so that there is no fictitious velocity at the next step. ---*/ - solver_container[selected_iZone][INST_0][MESH_0][MESH_SOL]->GetNodes()->Set_Solution_time_n(); - solver_container[selected_iZone][INST_0][MESH_0][MESH_SOL]->GetNodes()->Set_Solution_time_n1(); + solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->GetNodes()->Set_Solution_time_n(); + solver_container[selected_zone][INST_0][MESH_0][MESH_SOL]->GetNodes()->Set_Solution_time_n1(); } END_SU2_OMP_PARALLEL } @@ -143,7 +143,7 @@ void CDriver::BoundaryConditionsUpdate() { if (rank == MASTER_NODE) cout << "Updating boundary conditions." << endl; for (auto iZone = 0u; iZone < nZone; iZone++) { - geometry_container[iZone][INST_0][MESH_0]->UpdateCustomBoundaryConditions(geometry_container[iZone][INST_0],config_container[iZone]); + geometry_container[iZone][INST_0][MESH_0]->UpdateCustomBoundaryConditions(geometry_container[iZone][INST_0], config_container[iZone]); } } @@ -152,18 +152,14 @@ void CDriver::BoundaryConditionsUpdate() { //////////////////////////////////////////////////////////////////////////////// void CDriver::SetTranslationRate(passivedouble xDot, passivedouble yDot, passivedouble zDot) { - for (iZone = 0; iZone < nZone; iZone++) { - config_container[iZone]->SetTranslation_Rate(0, xDot); - config_container[iZone]->SetTranslation_Rate(1, yDot); - config_container[iZone]->SetTranslation_Rate(2, zDot); - } + main_config->SetTranslation_Rate(0, xDot); + main_config->SetTranslation_Rate(1, yDot); + main_config->SetTranslation_Rate(2, zDot); } void CDriver::SetRotationRate(passivedouble rot_x, passivedouble rot_y, passivedouble rot_z) { - for (iZone = 0; iZone < nZone; iZone++) { - config_container[iZone]->SetRotation_Rate(0, rot_x); - config_container[iZone]->SetRotation_Rate(1, rot_y); - config_container[iZone]->SetRotation_Rate(2, rot_z); - } + main_config->SetRotation_Rate(0, rot_x); + main_config->SetRotation_Rate(1, rot_y); + main_config->SetRotation_Rate(2, rot_z); } diff --git a/SU2_CFD/src/variables/CFEABoundVariable.cpp b/SU2_CFD/src/variables/CFEABoundVariable.cpp index 012498896ad..04d1e676260 100644 --- a/SU2_CFD/src/variables/CFEABoundVariable.cpp +++ b/SU2_CFD/src/variables/CFEABoundVariable.cpp @@ -65,8 +65,6 @@ void CFEABoundVariable::Set_FlowTraction_n() { FlowTraction_n = FlowTraction; } void CFEABoundVariable::Set_SurfaceLoad_Res_n() { Residual_Ext_Surf_n = Residual_Ext_Surf; } -void CFEABoundVariable::Clear_FlowTraction() { FlowTraction.setConstant(0.0); } - void CFEABoundVariable::Clear_SurfaceLoad_Res() { Residual_Ext_Surf.setConstant(0.0); } void CFEABoundVariable::RegisterFlowTraction(bool reset) { diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 81ed8e4f1c5..e86bf444a5f 100644 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -256,6 +256,8 @@ def main(): polar_naca0012.test_vals = [-1.217981, 4.256386, 0.009084, 0.016823] polar_naca0012.test_vals_aarch64 = [-1.718925, 3.711429, 0.009217, 0.007784] polar_naca0012.command = TestCase.Command(exec = "compute_polar.py", param = "-i 11") + # flaky test on arm64 + polar_naca0012.enabled_on_cpu_arch = ["x86_64"] test_list.append(polar_naca0012) # HYPERSONIC FLOW PAST BLUNT BODY @@ -1349,6 +1351,17 @@ def main(): pywrapper_fsi2d.multizone = True test_list.append(pywrapper_fsi2d) + # Unsteady FSI with custom load + pywrapper_unsteadyFSI = TestCase('pywrapper_unsteadyFSI') + pywrapper_unsteadyFSI.cfg_dir = "py_wrapper/dyn_fsi" + pywrapper_unsteadyFSI.cfg_file = "config.cfg" + pywrapper_unsteadyFSI.test_iter = 4 + pywrapper_unsteadyFSI.test_vals = [0, 31, 5, 58, -1.756780, -2.828276, -7.652558, -6.863929, 1.5618e-04] + pywrapper_unsteadyFSI.command = TestCase.Command("mpirun -np 2", "python", "run.py") + pywrapper_unsteadyFSI.unsteady = True + pywrapper_unsteadyFSI.multizone = True + test_list.append(pywrapper_unsteadyFSI) + # Unsteady CHT pywrapper_unsteadyCHT = TestCase('pywrapper_unsteadyCHT') pywrapper_unsteadyCHT.cfg_dir = "py_wrapper/flatPlate_unsteady_CHT" diff --git a/TestCases/py_wrapper/dyn_fsi/config.cfg b/TestCases/py_wrapper/dyn_fsi/config.cfg new file mode 100644 index 00000000000..28bafc8fccc --- /dev/null +++ b/TestCases/py_wrapper/dyn_fsi/config.cfg @@ -0,0 +1,31 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SU2 configuration file % +% Case description: Fluid Structure Interaction - Wall in channel % +% Author: R.Sanchez % +% Institution: Imperial College London % +% Date: 2016.02.01 % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SOLVER= MULTIPHYSICS + +CONFIG_LIST= (configFlow.cfg, configFEA.cfg) + +MARKER_ZONE_INTERFACE= (wallUpperF, wallUpperS, wallUpwF, wallUpwS, wallDownF, wallDownS) + +MULTIZONE_MESH= NO + +TIME_DOMAIN= YES +TIME_ITER= 5 +TIME_STEP= 0.002 +SCREEN_OUTPUT= (TIME_ITER, OUTER_ITER, INNER_ITER[0], INNER_ITER[1], DEFORM_ITER[0],\ + BGS_DENSITY[0], AVG_BGS_RES[1], RMS_DENSITY[0], RMS_UTOL[1], REFERENCE_NODE[1]) +RESTART_SOL= NO +RESTART_ITER= 2 + +OUTER_ITER= 10 +CONV_FIELD= BGS_DENSITY[0], AVG_BGS_RES[1] +CONV_RESIDUAL_MINVAL= -5 + +OUTPUT_FILES= RESTART, PARAVIEW +OUTPUT_WRT_FREQ= 1 +WRT_ZONE_HIST= YES diff --git a/TestCases/py_wrapper/dyn_fsi/configFEA.cfg b/TestCases/py_wrapper/dyn_fsi/configFEA.cfg new file mode 100644 index 00000000000..d782e3a1003 --- /dev/null +++ b/TestCases/py_wrapper/dyn_fsi/configFEA.cfg @@ -0,0 +1,54 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SU2 configuration file % +% Case description: Fluid Structure Interaction - Wall in channel % +% Author: R.Sanchez % +% Institution: Imperial College London % +% Date: 2016.02.01 % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SOLVER= ELASTICITY + +OBJECTIVE_FUNCTION= REFERENCE_NODE +REFERENCE_NODE= 3 +REFERENCE_NODE_DISPLACEMENT= (0.0, 0.0) +REFERENCE_NODE_PENALTY= 1.0 + +GEOMETRIC_CONDITIONS= LARGE_DEFORMATIONS +MATERIAL_MODEL= NEO_HOOKEAN +ELASTICITY_MODULUS= 1.5E5 +MATERIAL_DENSITY= 50 +FORMULATION_ELASTICITY_2D= PLANE_STRESS +POISSON_RATIO= 0.3 + +MARKER_CLAMPED= ( clamped ) +MARKER_FLUID_LOAD= ( wallDownS, wallUpperS, wallUpwS, internal ) + +BGS_RELAXATION= FIXED_PARAMETER +STAT_RELAX_PARAMETER= 0.5 + +LINEAR_SOLVER= CONJUGATE_GRADIENT +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-6 +LINEAR_SOLVER_ITER= 1000 + +TIME_DISCRE_FEA= NEWMARK_IMPLICIT +NEWMARK_BETA= 0.25 +NEWMARK_GAMMA= 0.5 + +PREDICTOR= YES +RELAXATION= YES + +INNER_ITER= 10 +CONV_FIELD= REL_RMS_UTOL +CONV_RESIDUAL_MINVAL= -5 + +SCREEN_WRT_FREQ_INNER= 100 + +MESH_FORMAT= SU2 +MESH_FILENAME= meshFEA.su2 + +VOLUME_FILENAME= results_wall +RESTART_FILENAME= restart_wall.dat +SOLUTION_FILENAME= solution_wall.dat + +HISTORY_OUTPUT= ITER, RMS_RES, STRUCT_COEFF, TAVG_STRUCT_COEFF diff --git a/TestCases/py_wrapper/dyn_fsi/configFlow.cfg b/TestCases/py_wrapper/dyn_fsi/configFlow.cfg new file mode 100644 index 00000000000..6b1df89cc7f --- /dev/null +++ b/TestCases/py_wrapper/dyn_fsi/configFlow.cfg @@ -0,0 +1,80 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SU2 configuration file % +% Case description: Fluid Structure Interaction - Wall in channel % +% Author: R.Sanchez % +% Institution: Imperial College London % +% Date: 2016.02.01 % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +SOLVER= NAVIER_STOKES +KIND_TURB_MODEL= NONE +DEFORM_MESH= YES + +MATH_PROBLEM= DIRECT + +AOA= 0.0 +REYNOLDS_NUMBER= 100 +REYNOLDS_LENGTH= 0.016 +MACH_NUMBER= 0.05 +MACH_MOTION= 0.2 + +INIT_OPTION= TD_CONDITIONS +FREESTREAM_OPTION= DENSITY_FS +FREESTREAM_DENSITY= 1.0 +FREESTREAM_PRESSURE= 17.85714286 +FREESTREAM_TEMPERATURE= 0.062207438 + +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.6E-2 + +REF_ORIGIN_MOMENT_X= 0 +REF_ORIGIN_MOMENT_Y= 0 +REF_ORIGIN_MOMENT_Z= 0 +REF_LENGTH= 0.016 +REF_AREA= 0.016 + +MARKER_INLET= ( inlet, 0.062207438, 23.36216288, 1.0, 0.0, 0.0 ) +MARKER_OUTLET= ( outlet, 17.85714286) +MARKER_EULER= ( upper, lower ) + +MARKER_HEATFLUX= ( wallUpwF, 0.0, wallDownF, 0.0, wallUpperF, 0.0 ) +MARKER_DEFORM_MESH= ( wallUpwF, wallDownF, wallUpperF ) +MARKER_DEFORM_MESH_SYM_PLANE= ( upper ) + +MARKER_PLOTTING= ( wallUpwF, wallDownF, wallUpperF ) +MARKER_MONITORING= ( wallUpwF, wallDownF, wallUpperF ) + +TIME_DISCRE_FLOW= EULER_IMPLICIT +CFL_NUMBER= 1000 +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 0.05 +LINEAR_SOLVER_ITER= 6 + +TIME_MARCHING= DUAL_TIME_STEPPING-2ND_ORDER +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES +SLOPE_LIMITER_FLOW= NONE + +DEFORM_LINEAR_SOLVER= CONJUGATE_GRADIENT +DEFORM_LINEAR_SOLVER_PREC= ILU +DEFORM_LINEAR_SOLVER_ERROR= 1E-9 +DEFORM_LINEAR_SOLVER_ITER= 100 +DEFORM_STIFFNESS_TYPE= WALL_DISTANCE + +INNER_ITER= 50 +CONV_FIELD= REL_RMS_DENSITY +CONV_RESIDUAL_MINVAL= -3 +CONV_STARTITER= 5 + +SCREEN_WRT_FREQ_INNER= 100 + +MESH_FORMAT= SU2 +MESH_FILENAME= meshFlow.su2 + +VOLUME_FILENAME= results_flow +RESTART_FILENAME= restart_flow.dat +SOLUTION_FILENAME= solution_flow.dat + +HISTORY_OUTPUT= ITER, RMS_RES, AERO_COEFF, TAVG_AERO_COEFF diff --git a/TestCases/py_wrapper/dyn_fsi/run.py b/TestCases/py_wrapper/dyn_fsi/run.py new file mode 100644 index 00000000000..ba6253db6b8 --- /dev/null +++ b/TestCases/py_wrapper/dyn_fsi/run.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python + +## \file run.py +# \brief Unsteady FSI case with custom load. +# \version 7.5.1 "Blackbird" +# +# SU2 Project Website: https://su2code.github.io +# +# The SU2 Project is maintained by the SU2 Foundation +# (http://su2foundation.org) +# +# Copyright 2012-2023, SU2 Contributors (cf. AUTHORS.md) +# +# SU2 is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License, or (at your option) any later version. +# +# SU2 is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with SU2. If not, see . + +import sys +import pysu2 +import math +from mpi4py import MPI + +def main(): + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + + # Initialize the corresponding driver of SU2, this includes solver preprocessing. + try: + driver = pysu2.CMultizoneDriver('config.cfg', 2, comm) + except TypeError as exception: + print('A TypeError occured in pysu2.CDriver : ', exception) + raise + + # By default the first zone is selected (flow in the case). + # Select the structural zone (second zone) to work with the FEA solver. + driver.SelectZone(1) + + # Get the ID of the marker where we will apply a force. + # This marker cannot be used by the fluid-solid interface otherwise the imposed + # load will be cleared when interpolating forces. + all_marker_ids = driver.GetMarkerIndices() + ctrl_id = all_marker_ids['internal'] if 'internal' in all_marker_ids else -1 + + # Number of vertices on the specified markers (per rank). + n_vertex_ctrl = driver.GetNumberMarkerNodes(ctrl_id) if ctrl_id >= 0 else 0 + + if rank == 0: + print("\n------------------------------ Begin Solver -----------------------------") + sys.stdout.flush() + + for time_iter in range(driver.GetNumberTimeIter()): + # Apply a custom load and then solve the time step. + time = time_iter * driver.GetUnsteadyTimeStep() + for i_vertex in range(n_vertex_ctrl): + i_point = driver.GetMarkerNode(ctrl_id, i_vertex) + if driver.GetNodeDomain(i_point): + driver.SetMarkerCustomFEALoad(ctrl_id, i_vertex, (-0.002 + 0.002 * math.cos(2 * math.pi * time / 0.02), 0)) + + driver.Preprocess(time_iter) + + driver.Run() + driver.Postprocess() + driver.Update() + + stop = driver.Monitor(time_iter) + driver.Output(time_iter) + if stop: break + + # Finalize the solver and exit cleanly. + driver.Finalize() + + +if __name__ == '__main__': + main() diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index d6512b98c8b..0ca25c6f3c4 100644 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -135,6 +135,8 @@ def main(): polar_naca0012.test_vals = [-1.243326, 4.224483, 0.016432, 0.016145] polar_naca0012.test_vals_aarch64 = [-1.811046, 3.612379, 0.012330, 0.009194] polar_naca0012.command = TestCase.Command(exec = "compute_polar.py", param = "-n 1 -i 11") + # flaky test on arm64 + polar_naca0012.enabled_on_cpu_arch = ["x86_64"] test_list.append(polar_naca0012) # HYPERSONIC FLOW PAST BLUNT BODY diff --git a/TestCases/vandv.py b/TestCases/vandv.py index baba0b84169..a5757cd7043 100644 --- a/TestCases/vandv.py +++ b/TestCases/vandv.py @@ -95,6 +95,7 @@ def main(): sandiajet_sst.cfg_file = "validation.cfg" sandiajet_sst.test_iter = 5 sandiajet_sst.test_vals = [-16.249917, -13.835991, -14.303372, -13.276035, -10.074262, -14.027223, 5, -1.672359, 5, -4.938477, 5, -3.462217, 2.5859e-04, 2.8215e-32, 4.5010e-68, 2.5859e-04, 4.0474e+03, 3.9468e+03, 4.9170e+01, 5.1441e+01] + sandiajet_sst.test_vals_aarch64 = [-16.249289, -13.833785, -14.303058, -13.276559, -10.267928, -14.027240, 5, -1.676412, 5, -4.815216, 5, -3.462247, 0.000259, 0, 0, 0.000259, 4047.4, 3946.8, 49.17, 51.441] test_list.append(sandiajet_sst) #################