diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 1860dea953f1..a79802b7fb6e 100644
--- a/Common/include/CConfig.hpp
+++ b/Common/include/CConfig.hpp
@@ -196,6 +196,7 @@ class CConfig {
nMarker_NearFieldBound, /*!< \brief Number of near field boundary markers. */
nMarker_ActDiskInlet, /*!< \brief Number of actuator disk inlet markers. */
nMarker_ActDiskOutlet, /*!< \brief Number of actuator disk outlet markers. */
+ nMarker_Deform_Mesh_Sym_Plane, /*!< \brief Number of markers with symmetric deformation */
nMarker_Deform_Mesh, /*!< \brief Number of deformable markers at the boundary. */
nMarker_Fluid_Load, /*!< \brief Number of markers in which the flow load is computed/employed. */
nMarker_Fluid_InterfaceBound, /*!< \brief Number of fluid interface markers. */
@@ -242,6 +243,7 @@ class CConfig {
*Marker_TurboBoundOut, /*!< \brief Turbomachinery performance boundary donor markers. */
*Marker_NearFieldBound, /*!< \brief Near Field boundaries markers. */
*Marker_Deform_Mesh, /*!< \brief Deformable markers at the boundary. */
+ *Marker_Deform_Mesh_Sym_Plane, /*!< \brief Marker with symmetric deformation. */
*Marker_Fluid_Load, /*!< \brief Markers in which the flow load is computed/employed. */
*Marker_Fluid_InterfaceBound, /*!< \brief Fluid interface markers. */
*Marker_CHTInterface, /*!< \brief Conjugate heat transfer interface markers. */
@@ -703,6 +705,7 @@ class CConfig {
*Marker_All_DV, /*!< \brief Global index for design variable markers using the grid information. */
*Marker_All_Moving, /*!< \brief Global index for moving surfaces using the grid information. */
*Marker_All_Deform_Mesh, /*!< \brief Global index for deformable markers at the boundary. */
+ *Marker_All_Deform_Mesh_Sym_Plane, /*!< \brief Global index for markers with symmetric deformations. */
*Marker_All_Fluid_Load, /*!< \brief Global index for markers in which the flow load is computed/employed. */
*Marker_All_PyCustom, /*!< \brief Global index for Python customizable surfaces using the grid information. */
*Marker_All_Designing, /*!< \brief Global index for moving using the grid information. */
@@ -717,6 +720,7 @@ class CConfig {
*Marker_CfgFile_MixingPlaneInterface, /*!< \brief Global index for MixingPlane interface using the config information. */
*Marker_CfgFile_Moving, /*!< \brief Global index for moving surfaces using the config information. */
*Marker_CfgFile_Deform_Mesh, /*!< \brief Global index for deformable markers at the boundary. */
+ *Marker_CfgFile_Deform_Mesh_Sym_Plane, /*!< \brief Global index for markers with symmetric deformations. */
*Marker_CfgFile_Fluid_Load, /*!< \brief Global index for markers in which the flow load is computed/employed. */
*Marker_CfgFile_PyCustom, /*!< \brief Global index for Python customizable surfaces using the config information. */
*Marker_CfgFile_DV, /*!< \brief Global index for design variable markers using the config information. */
@@ -3409,6 +3413,13 @@ class CConfig {
*/
void SetMarker_All_Deform_Mesh(unsigned short val_marker, unsigned short val_deform) { Marker_All_Deform_Mesh[val_marker] = val_deform; }
+ /*!
+ * \brief Set if a marker val_marker allows deformation at the boundary.
+ * \param[in] val_marker - Index of the marker in which we are interested.
+ * \param[in] val_interface - 0 or 1 depending if the the marker is or not a DEFORM_MESH_SYM_PLANE marker.
+ */
+ void SetMarker_All_Deform_Mesh_Sym_Plane(unsigned short val_marker, unsigned short val_deform) { Marker_All_Deform_Mesh_Sym_Plane[val_marker] = val_deform; }
+
/*!
* \brief Set if a in marker val_marker the flow load will be computed/employed.
* \param[in] val_marker - Index of the marker in which we are interested.
@@ -3546,6 +3557,13 @@ class CConfig {
*/
unsigned short GetMarker_All_Deform_Mesh(unsigned short val_marker) const { return Marker_All_Deform_Mesh[val_marker]; }
+ /*!
+ * \brief Get whether marker val_marker is a DEFORM_MESH_SYM_PLANE marker
+ * \param[in] val_marker - 0 or 1 depending if the the marker belongs to the DEFORM_MESH_SYM_PLANE subset.
+ * \return 0 or 1 depending if the marker belongs to the DEFORM_MESH_SYM_PLANE subset.
+ */
+ unsigned short GetMarker_All_Deform_Mesh_Sym_Plane(unsigned short val_marker) const { return Marker_All_Deform_Mesh_Sym_Plane[val_marker]; }
+
/*!
* \brief Get whether marker val_marker is a Fluid_Load marker
* \param[in] val_marker - 0 or 1 depending if the the marker belongs to the Fluid_Load subset.
@@ -6092,6 +6110,12 @@ class CConfig {
*/
unsigned short GetMarker_CfgFile_Deform_Mesh(string val_marker) const;
+ /*!
+ * \brief Get the DEFORM_MESH_SYM_PLANE information from the config definition for the marker val_marker.
+ * \return DEFORM_MESH_SYM_PLANE information of the boundary in the config information for the marker val_marker.
+ */
+ unsigned short GetMarker_CfgFile_Deform_Mesh_Sym_Plane(string val_marker) const;
+
/*!
* \brief Get the Fluid_Load information from the config definition for the marker val_marker.
* \return Fluid_Load information of the boundary in the config information for the marker val_marker.
@@ -6418,6 +6442,12 @@ class CConfig {
*/
unsigned short GetMarker_Deform_Mesh(string val_marker) const;
+ /*!
+ * \brief Get the internal index for a DEFORM_MESH_SYM_PLANE boundary val_marker.
+ * \return Internal index for a DEFORM_MESH_SYM_PLANE boundary val_marker.
+ */
+ unsigned short GetMarker_Deform_Mesh_Sym_Plane(string val_marker) const;
+
/*!
* \brief Get the internal index for a Fluid_Load boundary val_marker.
* \return Internal index for a Fluid_Load boundary val_marker.
@@ -6440,6 +6470,14 @@ class CConfig {
*/
string GetMarker_Deform_Mesh_TagBound(unsigned short val_marker) const { return Marker_Deform_Mesh[val_marker]; }
+ /*!
+ * \brief Get the name of the DEFORM_MESH_SYM_PLANE boundary defined in the geometry file.
+ * \param[in] val_marker - Value of the marker in which we are interested.
+ * \return Name that is in the geometry file for the surface that
+ * has the marker val_marker.
+ */
+ string GetMarker_Deform_Mesh_Sym_Plane_TagBound(unsigned short val_marker) const { return Marker_Deform_Mesh_Sym_Plane[val_marker]; }
+
/*!
* \brief Get the name of the Fluid_Load boundary defined in the geometry file.
* \param[in] val_marker - Value of the marker in which we are interested.
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index 782adcc18905..2d92146202af 100644
--- a/Common/src/CConfig.cpp
+++ b/Common/src/CConfig.cpp
@@ -789,6 +789,7 @@ void CConfig::SetPointersNull(void) {
Marker_CfgFile_MixingPlaneInterface = nullptr; Marker_All_MixingPlaneInterface = nullptr;
Marker_CfgFile_ZoneInterface = nullptr;
Marker_CfgFile_Deform_Mesh = nullptr; Marker_All_Deform_Mesh = nullptr;
+ Marker_CfgFile_Deform_Mesh_Sym_Plane = nullptr; Marker_All_Deform_Mesh_Sym_Plane = nullptr;
Marker_CfgFile_Fluid_Load = nullptr; Marker_All_Fluid_Load = nullptr;
Marker_CfgFile_Turbomachinery = nullptr; Marker_All_Turbomachinery = nullptr;
@@ -814,7 +815,7 @@ void CConfig::SetPointersNull(void) {
Marker_Euler = nullptr; Marker_FarField = nullptr; Marker_Custom = nullptr;
Marker_SymWall = nullptr; Marker_PerBound = nullptr;
Marker_PerDonor = nullptr; Marker_NearFieldBound = nullptr;
- Marker_Deform_Mesh = nullptr; Marker_Fluid_Load = nullptr;
+ Marker_Deform_Mesh = nullptr; Marker_Deform_Mesh_Sym_Plane= nullptr; Marker_Fluid_Load = nullptr;
Marker_Inlet = nullptr; Marker_Outlet = nullptr;
Marker_Supersonic_Inlet = nullptr; Marker_Supersonic_Outlet= nullptr; Marker_Smoluchowski_Maxwell = nullptr;
Marker_Isothermal = nullptr; Marker_HeatFlux = nullptr; Marker_EngineInflow = nullptr;
@@ -1387,6 +1388,8 @@ void CConfig::SetConfig_Options() {
addStringListOption("MARKER_FLUID_INTERFACE", nMarker_Fluid_InterfaceBound, Marker_Fluid_InterfaceBound);
/*!\brief MARKER_DEFORM_MESH\n DESCRIPTION: Deformable marker(s) at the interface \ingroup Config*/
addStringListOption("MARKER_DEFORM_MESH", nMarker_Deform_Mesh, Marker_Deform_Mesh);
+ /*!\brief MARKER_DEFORM_MESH_SYM_PLANE\n DESCRIPTION: Symmetry plane for mesh deformation only \ingroup Config*/
+ addStringListOption("MARKER_DEFORM_MESH_SYM_PLANE", nMarker_Deform_Mesh_Sym_Plane, Marker_Deform_Mesh_Sym_Plane);
/*!\brief MARKER_FLUID_LOAD\n DESCRIPTION: Marker(s) in which the flow load is computed/applied \ingroup Config*/
addStringListOption("MARKER_FLUID_LOAD", nMarker_Fluid_Load, Marker_Fluid_Load);
/*!\brief MARKER_FSI_INTERFACE \n DESCRIPTION: ZONE interface boundary marker(s) \ingroup Config*/
@@ -3707,7 +3710,7 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_
/*--- If it is not specified, set the mesh motion mach number
equal to the freestream value. ---*/
- if (GetGrid_Movement() && Mach_Motion == 0.0)
+ if (GetDynamic_Grid() && Mach_Motion == 0.0)
Mach_Motion = Mach;
/*--- Set the boolean flag if we are in a rotating frame (source term). ---*/
@@ -5042,7 +5045,7 @@ void CConfig::SetMarkers(unsigned short val_software) {
iMarker_Monitoring, iMarker_Designing, iMarker_GeoEval, iMarker_Plotting, iMarker_Analyze,
iMarker_DV, iMarker_Moving, iMarker_PyCustom, iMarker_Supersonic_Inlet, iMarker_Supersonic_Outlet,
iMarker_Clamped, iMarker_ZoneInterface, iMarker_CHTInterface, iMarker_Load_Dir, iMarker_Disp_Dir, iMarker_Load_Sine,
- iMarker_Fluid_Load, iMarker_Deform_Mesh,
+ iMarker_Fluid_Load, iMarker_Deform_Mesh, iMarker_Deform_Mesh_Sym_Plane,
iMarker_ActDiskInlet, iMarker_ActDiskOutlet,
iMarker_Turbomachinery, iMarker_MixingPlaneInterface;
@@ -5088,6 +5091,7 @@ void CConfig::SetMarkers(unsigned short val_software) {
Marker_All_DV = new unsigned short[nMarker_All] (); // Store whether the boundary should be affected by design variables.
Marker_All_Moving = new unsigned short[nMarker_All] (); // Store whether the boundary should be in motion.
Marker_All_Deform_Mesh = new unsigned short[nMarker_All] (); // Store whether the boundary is deformable.
+ Marker_All_Deform_Mesh_Sym_Plane = new unsigned short[nMarker_All] (); //Store wheter the boundary will follow the deformation
Marker_All_Fluid_Load = new unsigned short[nMarker_All] (); // Store whether the boundary computes/applies fluid loads.
Marker_All_PyCustom = new unsigned short[nMarker_All] (); // Store whether the boundary is Python customizable.
Marker_All_PerBound = new short[nMarker_All] (); // Store whether the boundary belongs to a periodic boundary.
@@ -5112,6 +5116,7 @@ void CConfig::SetMarkers(unsigned short val_software) {
Marker_CfgFile_DV = new unsigned short[nMarker_CfgFile] ();
Marker_CfgFile_Moving = new unsigned short[nMarker_CfgFile] ();
Marker_CfgFile_Deform_Mesh = new unsigned short[nMarker_CfgFile] ();
+ Marker_CfgFile_Deform_Mesh_Sym_Plane= new unsigned short[nMarker_CfgFile] ();
Marker_CfgFile_Fluid_Load = new unsigned short[nMarker_CfgFile] ();
Marker_CfgFile_PerBound = new unsigned short[nMarker_CfgFile] ();
Marker_CfgFile_Turbomachinery = new unsigned short[nMarker_CfgFile] ();
@@ -5511,6 +5516,13 @@ void CConfig::SetMarkers(unsigned short val_software) {
Marker_CfgFile_Deform_Mesh[iMarker_CfgFile] = YES;
}
+ for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) {
+ Marker_CfgFile_Deform_Mesh_Sym_Plane[iMarker_CfgFile] = NO;
+ for (iMarker_Deform_Mesh_Sym_Plane = 0; iMarker_Deform_Mesh_Sym_Plane < nMarker_Deform_Mesh_Sym_Plane; iMarker_Deform_Mesh_Sym_Plane++)
+ if (Marker_CfgFile_TagBound[iMarker_CfgFile] == Marker_Deform_Mesh_Sym_Plane[iMarker_Deform_Mesh_Sym_Plane])
+ Marker_CfgFile_Deform_Mesh_Sym_Plane[iMarker_CfgFile] = YES;
+ }
+
for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++) {
Marker_CfgFile_Fluid_Load[iMarker_CfgFile] = NO;
for (iMarker_Fluid_Load = 0; iMarker_Fluid_Load < nMarker_Fluid_Load; iMarker_Fluid_Load++)
@@ -5532,7 +5544,8 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) {
unsigned short iMarker_Euler, iMarker_Custom, iMarker_FarField,
iMarker_SymWall, iMarker_PerBound, iMarker_NearFieldBound,
iMarker_Fluid_InterfaceBound, iMarker_Inlet, iMarker_Riemann,
- iMarker_Deform_Mesh, iMarker_Fluid_Load, iMarker_Smoluchowski_Maxwell, iWall_Catalytic,
+ iMarker_Deform_Mesh, iMarker_Deform_Mesh_Sym_Plane, iMarker_Fluid_Load,
+ iMarker_Smoluchowski_Maxwell, iWall_Catalytic,
iMarker_Giles, iMarker_Outlet, iMarker_Isothermal, iMarker_HeatFlux,
iMarker_EngineInflow, iMarker_EngineExhaust, iMarker_Displacement, iMarker_Damper,
iMarker_Load, iMarker_FlowLoad, iMarker_Internal, iMarker_Monitoring,
@@ -6821,6 +6834,15 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) {
BoundaryTable.PrintFooter();
}
+ if (nMarker_Deform_Mesh_Sym_Plane != 0) {
+ BoundaryTable << "Symmetric deformable mesh boundary";
+ for (iMarker_Deform_Mesh_Sym_Plane = 0; iMarker_Deform_Mesh_Sym_Plane < nMarker_Deform_Mesh_Sym_Plane; iMarker_Deform_Mesh_Sym_Plane++) {
+ BoundaryTable << Marker_Deform_Mesh_Sym_Plane[iMarker_Deform_Mesh_Sym_Plane];
+ if (iMarker_Deform_Mesh_Sym_Plane < nMarker_Deform_Mesh_Sym_Plane-1) BoundaryTable << " ";
+ }
+ BoundaryTable.PrintFooter();
+ }
+
if (nMarker_Fluid_Load != 0) {
BoundaryTable << "Fluid loads boundary";
for (iMarker_Fluid_Load = 0; iMarker_Fluid_Load < nMarker_Fluid_Load; iMarker_Fluid_Load++) {
@@ -7342,6 +7364,13 @@ unsigned short CConfig::GetMarker_CfgFile_Deform_Mesh(string val_marker) const {
return Marker_CfgFile_Deform_Mesh[iMarker_CfgFile];
}
+unsigned short CConfig::GetMarker_CfgFile_Deform_Mesh_Sym_Plane(string val_marker) const {
+ unsigned short iMarker_CfgFile;
+ for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++)
+ if (Marker_CfgFile_TagBound[iMarker_CfgFile] == val_marker) break;
+ return Marker_CfgFile_Deform_Mesh_Sym_Plane[iMarker_CfgFile];
+}
+
unsigned short CConfig::GetMarker_CfgFile_Fluid_Load(string val_marker) const {
unsigned short iMarker_CfgFile;
for (iMarker_CfgFile = 0; iMarker_CfgFile < nMarker_CfgFile; iMarker_CfgFile++)
@@ -7516,6 +7545,9 @@ CConfig::~CConfig(void) {
delete[] Marker_CfgFile_Deform_Mesh;
delete[] Marker_All_Deform_Mesh;
+ delete[] Marker_CfgFile_Deform_Mesh_Sym_Plane;
+ delete[] Marker_All_Deform_Mesh_Sym_Plane;
+
delete[] Marker_CfgFile_Fluid_Load;
delete[] Marker_All_Fluid_Load;
@@ -7780,6 +7812,7 @@ CConfig::~CConfig(void) {
delete[] Marker_PerDonor;
delete[] Marker_NearFieldBound;
delete[] Marker_Deform_Mesh;
+ delete[] Marker_Deform_Mesh_Sym_Plane;
delete[] Marker_Fluid_Load;
delete[] Marker_Fluid_InterfaceBound;
delete[] Marker_Inlet;
@@ -8485,6 +8518,16 @@ unsigned short CConfig::GetMarker_Deform_Mesh(string val_marker) const {
return iMarker_Deform_Mesh;
}
+unsigned short CConfig::GetMarker_Deform_Mesh_Sym_Plane(string val_marker) const {
+ unsigned short iMarker_Deform_Mesh_Sym_Plane;
+
+ /*--- Find the marker for this interface boundary. ---*/
+ for (iMarker_Deform_Mesh_Sym_Plane = 0; iMarker_Deform_Mesh_Sym_Plane < nMarker_Deform_Mesh_Sym_Plane; iMarker_Deform_Mesh_Sym_Plane++)
+ if (Marker_Deform_Mesh_Sym_Plane[iMarker_Deform_Mesh_Sym_Plane] == val_marker) break;
+
+ return iMarker_Deform_Mesh_Sym_Plane;
+}
+
unsigned short CConfig::GetMarker_Fluid_Load(string val_marker) const {
unsigned short iMarker_Fluid_Load;
diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp
index ac8e34710649..1a5b36a1f3ab 100644
--- a/Common/src/geometry/CPhysicalGeometry.cpp
+++ b/Common/src/geometry/CPhysicalGeometry.cpp
@@ -3593,6 +3593,7 @@ void CPhysicalGeometry::SetBoundaries(CConfig *config) {
config->SetMarker_All_DV(iMarker, config->GetMarker_CfgFile_DV(Marker_Tag));
config->SetMarker_All_Moving(iMarker, config->GetMarker_CfgFile_Moving(Marker_Tag));
config->SetMarker_All_Deform_Mesh(iMarker, config->GetMarker_CfgFile_Deform_Mesh(Marker_Tag));
+ config->SetMarker_All_Deform_Mesh_Sym_Plane(iMarker, config->GetMarker_CfgFile_Deform_Mesh_Sym_Plane(Marker_Tag));
config->SetMarker_All_Fluid_Load(iMarker, config->GetMarker_CfgFile_Fluid_Load(Marker_Tag));
config->SetMarker_All_PyCustom(iMarker, config->GetMarker_CfgFile_PyCustom(Marker_Tag));
config->SetMarker_All_PerBound(iMarker, config->GetMarker_CfgFile_PerBound(Marker_Tag));
@@ -3615,6 +3616,7 @@ void CPhysicalGeometry::SetBoundaries(CConfig *config) {
config->SetMarker_All_DV(iMarker, NO);
config->SetMarker_All_Moving(iMarker, NO);
config->SetMarker_All_Deform_Mesh(iMarker, NO);
+ config->SetMarker_All_Deform_Mesh_Sym_Plane(iMarker, NO);
config->SetMarker_All_Fluid_Load(iMarker, NO);
config->SetMarker_All_PyCustom(iMarker, NO);
config->SetMarker_All_PerBound(iMarker, NO);
@@ -4016,6 +4018,7 @@ void CPhysicalGeometry::LoadUnpartitionedSurfaceElements(CConfig *config,
config->SetMarker_All_DV(iMarker, config->GetMarker_CfgFile_DV(Marker_Tag));
config->SetMarker_All_Moving(iMarker, config->GetMarker_CfgFile_Moving(Marker_Tag));
config->SetMarker_All_Deform_Mesh(iMarker, config->GetMarker_CfgFile_Deform_Mesh(Marker_Tag));
+ config->SetMarker_All_Deform_Mesh_Sym_Plane(iMarker, config->GetMarker_CfgFile_Deform_Mesh_Sym_Plane(Marker_Tag));
config->SetMarker_All_Fluid_Load(iMarker, config->GetMarker_CfgFile_Fluid_Load(Marker_Tag));
config->SetMarker_All_PyCustom(iMarker, config->GetMarker_CfgFile_PyCustom(Marker_Tag));
config->SetMarker_All_PerBound(iMarker, config->GetMarker_CfgFile_PerBound(Marker_Tag));
diff --git a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp
index 635e5739134a..055a35527dfa 100644
--- a/SU2_CFD/include/drivers/CSinglezoneDriver.hpp
+++ b/SU2_CFD/include/drivers/CSinglezoneDriver.hpp
@@ -7,7 +7,7 @@
*
* SU2 Project Website: https://su2code.github.io
*
- * The SU2 Project is maintained by the SU2 Foundation
+ * The SU2 Project is maintained by the SU2 Foundation
* (http://su2foundation.org)
*
* Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
@@ -92,6 +92,11 @@ class CSinglezoneDriver : public CDriver {
*/
void DynamicMeshUpdate(unsigned long TimeIter) override;
+ /*!
+ * \brief Perform a mesh deformation as initial condition.
+ */
+ void SetInitialMesh() override;
+
/*!
* \brief Monitor
* \param ExtIter
diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp
index 96494dc069e9..f783e4938355 100644
--- a/SU2_CFD/src/iteration/CFluidIteration.cpp
+++ b/SU2_CFD/src/iteration/CFluidIteration.cpp
@@ -235,7 +235,6 @@ void CFluidIteration::Postprocess(COutput* output, CIntegration**** integration,
CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement,
CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) {
/*--- Temporary: enable only for single-zone driver. This should be removed eventually when generalized. ---*/
-
if (config[val_iZone]->GetSinglezone_Driver()) {
/*--- Compute the tractions at the vertices ---*/
solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->ComputeVertexTractions(geometry[val_iZone][val_iInst][MESH_0],
diff --git a/SU2_CFD/src/output/COutput.cpp b/SU2_CFD/src/output/COutput.cpp
index d3606ac58ed8..f56bab1015ae 100644
--- a/SU2_CFD/src/output/COutput.cpp
+++ b/SU2_CFD/src/output/COutput.cpp
@@ -1992,7 +1992,9 @@ bool COutput::WriteHistoryFile_Output(CConfig *config) {
}
bool COutput::WriteVolume_Output(CConfig *config, unsigned long Iter, bool force_writing){
- if (config->GetTime_Domain()) return ((Iter % config->GetVolume_Wrt_Freq() == 0)) || force_writing;
+ if (config->GetTime_Domain()){
+ return ((Iter % config->GetVolume_Wrt_Freq() == 0)) || force_writing;
+ }
else {
return ((Iter > 0) && (Iter % config->GetVolume_Wrt_Freq() == 0)) || force_writing;
}
diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp
index ec823ff9d19b..dfa52ad48e4c 100644
--- a/SU2_CFD/src/python_wrapper_structure.cpp
+++ b/SU2_CFD/src/python_wrapper_structure.cpp
@@ -27,6 +27,7 @@
#include "../include/drivers/CDriver.hpp"
+ #include "../include/drivers/CSinglezoneDriver.hpp"
void CDriver::PythonInterface_Preprocessing(CConfig **config, CGeometry ****geometry, CSolver *****solver){
@@ -891,6 +892,32 @@ void CFluidDriver::SetInitialMesh() {
//}
}
+void CSinglezoneDriver::SetInitialMesh() {
+
+ DynamicMeshUpdate(0);
+
+ SU2_OMP_PARALLEL {
+ // Overwrite fictious velocities
+ for (iMesh = 0u; iMesh <= config_container[ZONE_0]->GetnMGLevels(); iMesh++) {
+ SU2_OMP_FOR_STAT(roundUpDiv(geometry_container[ZONE_0][INST_0][iMesh]->GetnPoint(),omp_get_max_threads()))
+ for (unsigned long iPoint = 0; iPoint < geometry_container[ZONE_0][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[ZONE_0][INST_0][iMesh]->nodes->SetGridVel(iPoint, Grid_Vel);
+ }
+ /*--- Push back the volume. ---*/
+ geometry_container[ZONE_0][INST_0][iMesh]->nodes->SetVolume_n();
+ geometry_container[ZONE_0][INST_0][iMesh]->nodes->SetVolume_nM1();
+ }
+ /*--- Push back the solution so that there is no fictious velocity at the next step. ---*/
+ solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->Set_Solution_time_n();
+ solver_container[ZONE_0][INST_0][MESH_0][MESH_SOL]->GetNodes()->Set_Solution_time_n1();
+ }
+}
+
void CFluidDriver::SetVertexTtotal(unsigned short iMarker, unsigned long iVertex, passivedouble val_Ttotal_passive){
su2double val_Ttotal = val_Ttotal_passive;
@@ -1231,4 +1258,3 @@ void CDriver::SetInlet_Angle(unsigned short iMarker, passivedouble alpha){
}
}
-
diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp
index bc35084d91eb..8fce143bfd0d 100644
--- a/SU2_CFD/src/solvers/CFEASolver.cpp
+++ b/SU2_CFD/src/solvers/CFEASolver.cpp
@@ -1788,11 +1788,9 @@ void CFEASolver::BC_Clamped_Post(CGeometry *geometry, CNumerics *numerics, const
void CFEASolver::BC_Sym_Plane(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker) {
if (geometry->GetnElem_Bound(val_marker) == 0) return;
-
const bool dynamic = config->GetTime_Domain();
/*--- Determine axis of symmetry based on the normal of the first element in the marker. ---*/
-
const su2double* nodeCoord[MAXNNODE_2D] = {nullptr};
const bool quad = (geometry->bound[val_marker][0]->GetVTK_Type() == QUADRILATERAL);
@@ -1802,7 +1800,6 @@ void CFEASolver::BC_Sym_Plane(CGeometry *geometry, CNumerics *numerics, const CC
auto iPoint = geometry->bound[val_marker][0]->GetNode(iNode);
nodeCoord[iNode] = geometry->nodes->GetCoord(iPoint);
}
-
su2double normal[MAXNDIM] = {0.0};
switch (nNodes) {
@@ -1828,7 +1825,6 @@ void CFEASolver::BC_Sym_Plane(CGeometry *geometry, CNumerics *numerics, const CC
/*--- Set and enforce solution at current and previous time-step ---*/
nodes->SetSolution(iPoint, axis, 0.0);
-
if (dynamic) {
nodes->SetSolution_Vel(iPoint, axis, 0.0);
nodes->SetSolution_Accel(iPoint, axis, 0.0);
@@ -1839,9 +1835,10 @@ void CFEASolver::BC_Sym_Plane(CGeometry *geometry, CNumerics *numerics, const CC
/*--- Set and enforce 0 solution for mesh deformation ---*/
nodes->SetBound_Disp(iPoint, axis, 0.0);
-
LinSysSol(iPoint, axis) = 0.0;
- LinSysReact(iPoint, axis) = 0.0;
+ if (LinSysReact.GetLocSize() > 0){
+ LinSysReact(iPoint, axis) = 0.0;
+ }
Jacobian.EnforceSolutionAtDOF(iPoint, axis, su2double(0.0), LinSysRes);
}
@@ -2708,7 +2705,6 @@ void CFEASolver::GeneralizedAlpha_UpdateLoads(CGeometry *geometry, const CConfig
void CFEASolver::Solve_System(CGeometry *geometry, CConfig *config) {
/*--- Enforce solution at some halo points possibly not covered by essential BC markers. ---*/
-
Jacobian.InitiateComms(LinSysSol, geometry, config, SOLUTION_MATRIX);
Jacobian.CompleteComms(LinSysSol, geometry, config, SOLUTION_MATRIX);
@@ -2725,6 +2721,7 @@ void CFEASolver::Solve_System(CGeometry *geometry, CConfig *config) {
/*--- Solve or smooth the linear system. ---*/
auto iter = System.Solve(Jacobian, LinSysRes, LinSysSol, geometry, config);
+
SU2_OMP_MASTER
{
SetIterLinSolver(iter);
diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp
index f964cd51cd87..c124fa66cfa2 100644
--- a/SU2_CFD/src/solvers/CMeshSolver.cpp
+++ b/SU2_CFD/src/solvers/CMeshSolver.cpp
@@ -649,6 +649,7 @@ void CMeshSolver::SetBoundaryDisplacements(CGeometry *geometry, CNumerics *numer
/*--- Exceptions: symmetry plane, the receive boundaries and periodic boundaries should get a different treatment. ---*/
for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
if ((config->GetMarker_All_Deform_Mesh(iMarker) == NO) &&
+ (config->GetMarker_All_Deform_Mesh_Sym_Plane(iMarker) == NO) &&
(config->GetMarker_All_Moving(iMarker) == NO) &&
(config->GetMarker_All_KindBC(iMarker) != SYMMETRY_PLANE) &&
(config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) &&
@@ -661,6 +662,7 @@ void CMeshSolver::SetBoundaryDisplacements(CGeometry *geometry, CNumerics *numer
/*--- Symmetry plane is clamped, for now. ---*/
for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
if ((config->GetMarker_All_Deform_Mesh(iMarker) == NO) &&
+ (config->GetMarker_All_Deform_Mesh_Sym_Plane(iMarker) == NO) &&
(config->GetMarker_All_Moving(iMarker) == NO) &&
(config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE)) {
@@ -668,6 +670,7 @@ void CMeshSolver::SetBoundaryDisplacements(CGeometry *geometry, CNumerics *numer
}
}
+
/*--- Impose displacement boundary conditions. ---*/
for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
if ((config->GetMarker_All_Deform_Mesh(iMarker) == YES) ||
@@ -677,6 +680,15 @@ void CMeshSolver::SetBoundaryDisplacements(CGeometry *geometry, CNumerics *numer
}
}
+
+ /*--- Symmetry deform plane is not clamped ---*/
+ for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+ if (config->GetMarker_All_Deform_Mesh_Sym_Plane(iMarker) == YES) {
+
+ BC_Sym_Plane(geometry, numerics, config, iMarker);
+ }
+ }
+
/*--- Clamp far away nodes according to deform limit. ---*/
if ((config->GetDeform_Stiffness_Type() == SOLID_WALL_DISTANCE) &&
(config->GetDeform_Limit() < MaxDistance)) {
diff --git a/SU2_PY/FSI/PitchPlungeAirfoilStructuralTester.py b/SU2_PY/FSI/PitchPlungeAirfoilStructuralTester.py
deleted file mode 100755
index c670d5b513a1..000000000000
--- a/SU2_PY/FSI/PitchPlungeAirfoilStructuralTester.py
+++ /dev/null
@@ -1,712 +0,0 @@
-#!/usr/bin/env python
-
-## \file SolidSolverTester.py
-# \brief Structural solver tester (one or two degree of freedom) used for testing the Py wrapper for external FSI coupling.
-# \author David Thomas
-# \version 7.0.8 "Blackbird"
-#
-# SU2 Project Website: https://su2code.github.io
-#
-# The SU2 Project is maintained by the SU2 Foundation
-# (http://su2foundation.org)
-#
-# Copyright 2012-2020, 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 .
-
-# ----------------------------------------------------------------------
-# Imports
-# ----------------------------------------------------------------------
-
-import os, sys, shutil, copy
-import numpy as np
-import scipy as sp
-import scipy.linalg as linalg
-from math import *
-from util import switch
-
-# ----------------------------------------------------------------------
-# Config class
-# ----------------------------------------------------------------------
-
-class Point:
- """ Description. """
-
- def __init__(self):
- self.Coord0 = np.zeros((3,1))
- self.Coord = np.zeros((3,1))
- self.Coord_n = np.zeros((3,1))
- self.Vel = np.zeros((3,1))
- self.Vel_n = np.zeros((3,1))
- self.Force = np.zeros((3,1))
-
- def GetCoord0(self):
- return self.Coord0
-
- def GetCoord(self):
- return self.Coord
-
- def GetCoord_n(self):
- return self.Coord_n
-
- def GetVel(self):
- return self.Vel
-
- def GetVel_n(self):
- return self.Vel_n
-
- def GetForce(self):
- return self.Force
-
- def SetCoord0(self, val_Coord):
- x, y, z = val_Coord
- self.Coord0[0] = x
- self.Coord0[1] = y
- self.Coord0[2] = z
-
- def SetCoord(self, val_Coord):
- x, y, z = val_Coord
- self.Coord[0] = x
- self.Coord[1] = y
- self.Coord[2] = z
-
- def SetCoord_n(self, val_Coord):
- x, y, z = val_Coord
- self.Coord_n[0] = x
- self.Coord_n[1] = y
- self.Coord_n[2] = z
-
- def SetVel(self, val_Vel):
- vx, vy, vz = val_Vel
- self.Vel[0] = vx
- self.Vel[1] = vy
- self.Vel[2] = vz
-
- def SetVel_n(self, val_Vel):
- vx, vy, vz = val_Vel
- self.Vel_n[0] = vx
- self.Vel_n[1] = vy
- self.Vel_n[2] = vz
-
- def SetForce(self, val_Force):
- fx, fy, fz = val_Force
- self.Force[0] = fx
- self.Force[1] = fy
- self.Force[2] = fz
-
- def updateCoordVel(self):
- self.Coord_n = np.copy(self.Coord)
- self.Vel_n = np.copy(self.Vel)
-
-class Solver:
- """Description"""
-
- def __init__(self, config_fileName):
- """ Description. """
-
- self.Config_file = config_fileName
- self.Config = {}
-
- print("\n------------------------------ Configuring the structural tester solver for FSI simulation ------------------------------")
- self.__readConfig()
-
- self.Mesh_file = self.Config['MESH_FILE']
- self.FSI_marker = self.Config['MOVING_MARKER']
- self.Unsteady = (self.Config['TIME_MARCHING']=="YES")
- if self.Unsteady:
- print('Dynamic computation.')
- if self.Config['STRUCT_TYPE'] == "AIRFOIL":
- self.nDof = 2
- print("Structural model : pitching-plunging airfoil.")
- else:
- self.nDof = 0
-
-
- # Structural properties
- self.m = self.Config['SPRING_MASS'] # airfoil mass [kg]
- self.Kh = self.Config['SPRING_STIFFNESS'] # plunging stiffness [N/m]
- self.Ka = self.Config['TORSIONAL_STIFFNESS'] # pitching stiffness [N]
- self.If = self.Config['INERTIA_FLEXURAL'] # inertia around flexural axis [kg m^2]
- self.Ch = self.Config['SPRING_DAMPING'] # plunging damping [Ns/m]
- self.Ca = self.Config['TORSIONAL_DAMPING'] # pitching damping [Ns]
- self.c = self.Config['CORD'] # airfoil cord [m]
- self.b = self.c/2.0 # airfoil semi cord [m]
- self.xf = self.Config['FLEXURAL_AXIS'] # position of the flexural axis [m]
- self.xCG = self.Config['GRAVITY_CENTER'] # position of the center of gravity [m]
- self.S = self.m*(self.xCG - self.xf)
- self.h0 = self.Config['INITIAL_DISP']
- self.a0 = self.Config['INITIAL_ANGLE']
- self.startTime = self.Config['START_TIME']
- self.stopTime = self.Config['STOP_TIME']
- self.deltaT = self.Config['DELTA_T']
- self.rhoAlphaGen = self.Config['RHO']
-
- self.nDim= int()
- self.nElem = int()
- self.nPoint = int()
- self.nMarker = int()
- self.node = []
- self.markers = {}
-
- print("\n------------------------------ Reading the SU2 mesh ------------------------------")
- self.__readSU2Mesh()
-
- print("\n------------------------------ Creating the structural model ------------------------------")
- self.__setStructuralMatrices()
-
- print("\n------------------------------ Setting the integration parameters ------------------------------")
- self.__setIntegrationParameters()
- self.__setInitialConditions()
-
- def __readConfig(self):
- """ Description. """
-
- with open(self.Config_file) as configfile:
- while 1:
- line = configfile.readline()
- if not line:
- break
-
- # remove line returns
- line = line.strip('\r\n')
- # make sure it has useful data
- if (not "=" in line) or (line[0] == '%'):
- continue
- # split across equal sign
- line = line.split("=",1)
- this_param = line[0].strip()
- this_value = line[1].strip()
-
- for case in switch(this_param):
- #integer values
- #if case("NB_FSI_ITER") :
- #self.Config[this_param] = int(this_value)
- #break
-
- #float values
- if case("DELTA_T") : pass
- if case("START_TIME") : pass
- if case("STOP_TIME") : pass
- if case("SPRING_MASS") : pass
- if case("INERTIA_FLEXURAL") : pass
- if case("SPRING_STIFFNESS") : pass
- if case("SPRING_DAMPING") : pass
- if case("TORSIONAL_STIFFNESS") : pass
- if case("TORSIONAL_DAMPING") : pass
- if case("CORD") : pass
- if case("FLEXURAL_AXIS") : pass
- if case("GRAVITY_CENTER") : pass
- if case("INITIAL_DISP") : pass
- if case("INITIAL_ANGLE") : pass
- if case("RHO") :
- self.Config[this_param] = float(this_value)
- break
-
- #string values
- if case("TIME_MARCHING") : pass
- if case("MESH_FILE") : pass
- if case("CSD_SOLVER") : pass
- if case("MOVING_MARKER") : pass
- if case("STRUCT_TYPE") :
- self.Config[this_param] = this_value
- break
-
- if case():
- print(this_param + " is an invalid option !")
- break
-
- def __readSU2Mesh(self):
- """ Description. """
-
- with open(self.Mesh_file, 'r') as meshfile:
- print('Opened mesh file ' + self.Mesh_file + '.')
- while 1:
- line = meshfile.readline()
- if not line:
- break
-
- pos = line.find('NDIM')
- if pos != -1:
- line = line.strip('\r\n')
- line = line.split("=",1)
- self.nDim = int(line[1])
- continue
-
- pos = line.find('NELEM')
- if pos != -1:
- line = line.strip('\r\n')
- line = line.split("=",1)
- self.nElem = int(line[1])
- continue
-
- pos = line.find('NPOIN')
- if pos != -1:
- line = line.strip('\r\n')
- line = line.split("=",1)
- self.nPoint = int(line[1])
- for iPoint in range(self.nPoint):
- self.node.append(Point())
- line = meshfile.readline()
- line = line.strip('\r\n')
- line = line.split(' ',self.nDim)
- x = float(line[0])
- y = float(line[1])
- z = 0.0
- if self.nDim == 3:
- z = float(line[2])
- self.node[iPoint].SetCoord((x,y,z))
- self.node[iPoint].SetCoord0((x,y,z))
- self.node[iPoint].SetCoord_n((x,y,z))
- continue
-
- pos = line.find('NMARK')
- if pos != -1:
- line = line.strip('\r\n')
- line = line.split("=",1)
- self.nMarker = int(line[1])
- continue
-
- pos = line.find('MARKER_TAG')
- if pos != -1:
- line = line.strip('\r\n')
- line = line.replace(" ", "")
- line = line.split("=",1)
- markerTag = line[1]
- if markerTag == self.FSI_marker:
- self.markers[markerTag] = []
- line = meshfile.readline()
- line = line.strip('\r\n')
- line = line.split("=",1)
- nElem = int(line[1])
- for iElem in range(nElem):
- line = meshfile.readline()
- line = line.strip('\r\n')
- line = line.split(' ',1)
- elemType = int(line[0])
- if elemType == 3:
- nodes = line[1].split(' ', 1)
- if not int(nodes[0]) in self.markers[markerTag]:
- self.markers[markerTag].append(int(nodes[0]))
- if not int(nodes[1]) in self.markers[markerTag]:
- self.markers[markerTag].append(int(nodes[1]))
- else:
- print("Element type {} is not recognized !!".format(elemType))
- continue
- else:
- continue
-
- print("Number of dimensions: {}".format(self.nDim))
- print("Number of elements: {}".format(self.nElem))
- print("Number of point: {}".format(self.nPoint))
- print("Number of markers: {}".format(self.nMarker))
- if len(self.markers) > 0:
- print("Moving marker(s):")
- for mark in self.markers.keys():
- print(mark)
-
- def __setStructuralMatrices(self):
- """ Descriptions. """
-
- self.M = np.zeros((self.nDof, self.nDof))
- self.K = np.zeros((self.nDof, self.nDof))
- self.C = np.zeros((self.nDof, self.nDof))
-
- self.q = np.zeros((self.nDof, 1))
- self.qdot = np.zeros((self.nDof, 1))
- self.qddot = np.zeros((self.nDof, 1))
- self.a = np.zeros((self.nDof, 1))
-
- self.q_n = np.zeros((self.nDof, 1))
- self.qdot_n = np.zeros((self.nDof, 1))
- self.qddot_n = np.zeros((self.nDof, 1))
- self.a_n = np.zeros((self.nDof, 1))
-
- self.F = np.zeros((self.nDof, 1))
-
- if self.Config['STRUCT_TYPE'] == "AIRFOIL":
- print('Setting pitching-plunging airfoil system')
- print('Number of DOF : ')
-
- self.centerOfRotation = np.zeros((3,1))
- self.centerOfRotation[0] = self.xf
- self.centerOfRotation_n = np.zeros((3,1))
- self.centerOfRotation_n[0] = self.xf
- self.M[0][0] = self.m
- self.M[0][1] = self.S
- self.M[1][0] = self.S
- self.M[1][1] = self.If
- self.C[0][0] = self.Ch
- self.C[1][1] = self.Ca
- self.K[0][0] = self.Kh
- self.K[1][1] = self.Ka
-
- print('Airfoil mass : {} [kg]'.format(self.m))
- print('Airfoil cord : {} [m]'.format(self.c))
- print('Position of the flexural axis (from the leading edge) : {} [m]'.format(self.xf))
- print('Position of the center of gravity (from the leading edge) : {} [m]'.format(self.xCG))
- print('Inertia around the flexural axis : {} [kg m^2]'.format(self.If))
- print('Static unbalance : {} [kg m]'.format(self.S))
- print('Plunging stiffness : {} [N/m]'.format(self.Kh))
- print('Plunging dampimg : {} [Ns/m]'.format(self.Ch))
- print('Pitching stiffness : {} [N]'.format(self.Ka))
- print('Pitching dampimg : {} [Ns]'.format(self.Ca))
-
- def __setIntegrationParameters(self):
- """ Description. """
-
- self.alpha_m = (2.0*self.rhoAlphaGen-1.0)/(self.rhoAlphaGen+1.0)
- self.alpha_f = (self.rhoAlphaGen)/(self.rhoAlphaGen+1.0)
- self.gamma = 0.5+self.alpha_f-self.alpha_m
- self.beta = 0.25*(self.gamma+0.5)**2
-
- self.gammaPrime = self.gamma/(self.deltaT*self.beta)
- self.betaPrime = (1.0-self.alpha_m)/((self.deltaT**2)*self.beta*(1.0-self.alpha_f))
-
- print('Time integration with the alpha-generalized algorithm.')
- print('rho : {}'.format(self.rhoAlphaGen))
- print('alpha_m : {}'.format(self.alpha_m))
- print('alpha_f : {}'.format(self.alpha_f))
- print('gamma : {}'.format(self.gamma))
- print('beta : {}'.format(self.beta))
- print('gammaPrime : {}'.format(self.gammaPrime))
- print('betaPrime : {}'.format(self.betaPrime))
-
- def __setInitialConditions(self):
- """ Description. """
-
- print('Setting initial conditions.')
- self.__reset(self.F)
- self.__reset(self.q_n)
-
- print('Initial plunge displacement : {} [m]'.format(self.h0))
- self.q[0] = self.h0
- if self.nDof == 2:
- print('Initial pitch angle : {} [rad]'.format(self.a0))
- self.q[1] = self.a0
-
- RHS = np.zeros((self.nDof,1))
- RHS += self.F
- RHS -= self.C.dot(self.qdot)
- RHS -= self.K.dot(self.q)
- self.qddot = linalg.solve(self.M, RHS)
- self.a = np.copy(self.qddot)
-
- self.centerOfRotation[1] = self.q[0]
-
- def __reset(self, vector):
- """ Description. """
-
- for ii in range(vector.shape[0]):
- vector[ii] = 0.0
-
- def __computeInterfacePosVel(self, initialize):
- """ Description. """
-
- dTheta = 0.0
- dPhi = 0.0
- newCenter = np.zeros((3,1))
- Centerdot = np.zeros((3,1))
- newVel = np.zeros((3,1))
-
- if self.nDof == 2:
- dPsi = -(self.q[1] - self.q_n[1])
- newCenter[0] = self.centerOfRotation[0]
- newCenter[1] = -self.q[0]
- newCenter[2] = self.centerOfRotation[2]
- Centerdot[0] = 0.0
- Centerdot[1] = -self.qdot[0]
- Centerdot[2] = 0.0
- psidot = self.qdot[1]
-
- cosPsi = cos(dPsi)
- sinPsi = sin(dPsi)
-
- rotMatrix = np.array([[cosPsi, -sinPsi, 0.0],[sinPsi, cosPsi, 0.0],[0.0, 0.0, 1.0]])
-
- for iMarker in self.markers.keys():
- vertexList = self.markers[iMarker]
- for iPoint in vertexList:
- Coord = self.node[iPoint].GetCoord()
- Coord_n = self.node[iPoint].GetCoord_n()
-
- if self.Unsteady:
- r = Coord_n - self.centerOfRotation_n
- else:
- r = Coord - self.centerOfRotation
-
- rotCoord = rotMatrix.dot(r)
-
- newCoord = newCenter + rotCoord
- newVel[0] = Centerdot[0]+psidot*(newCoord[1]-newCenter[1])
- newVel[1] = Centerdot[1]-psidot*(newCoord[0]-newCenter[0])
- newVel[2] = Centerdot[2]+0.0
-
- self.node[iPoint].SetCoord((newCoord[0], newCoord[1], newCoord[2]))
- self.node[iPoint].SetVel((newVel[0], newVel[1], newVel[2]))
-
- if initialize:
- self.node[iPoint].SetCoord_n((newCoord[0], newCoord[1], newCoord[2]))
- self.node[iPoint].SetVel_n((newVel[0], newVel[1], newVel[2]))
-
- self.centerOfRotation = np.copy(newCenter)
-
- def __temporalIteration(self):
- """ Description. """
-
- eps = 1e-6
-
- self.__SetLoads()
-
- # Prediction step
- self.__reset(self.qddot)
- self.__reset(self.a)
-
- self.a += (self.alpha_f)/(1-self.alpha_m)*self.qddot_n
- self.a -= (self.alpha_m)/(1-self.alpha_m)*self.a_n
-
- self.q = np.copy(self.q_n)
- self.q += self.deltaT*self.qdot_n
- self.q += (0.5-self.beta)*self.deltaT*self.deltaT*self.a_n
- self.q += self.deltaT*self.deltaT*self.beta*self.a
-
- self.qdot = np.copy(self.qdot_n)
- self.qdot += (1-self.gamma)*self.deltaT*self.a_n
- self.qdot += self.deltaT*self.gamma*self.a
-
- # Correction step
- res = self.__ComputeResidual()
-
- while linalg.norm(res) >= eps:
- St = self.__TangentOperator()
- Deltaq = -1*(linalg.solve(St,res))
- self.q += Deltaq
- self.qdot += self.gammaPrime*Deltaq
- self.qddot += self.betaPrime*Deltaq
- res = self.__ComputeResidual()
-
- self.a += (1-self.alpha_f)/(1-self.alpha_m)*self.qddot
-
-
- def __SetLoads(self):
- """ Description """
-
- makerID = self.markers.keys()[0]
- nodeList = self.markers[makerID]
-
- FX = 0.0
- FY = 0.0
- FZ = 0.0
- MZ = 0.0
-
- for iPoint in nodeList:
- Force = self.node[iPoint].GetForce()
- Coord = self.node[iPoint].GetCoord()
- FX += float(Force[0])
- FY += float(Force[1])
- FZ += float(Force[2])
- MZ += float(Force[1]*(Coord[0]-self.centerOfRotation[0])-Force[0]*(Coord[1]-self.centerOfRotation[1]))
-
- self.F[0] = -FY
- self.F[1] = -MZ
-
- def __ComputeResidual(self):
- """ Description. """
-
- res = self.M.dot(self.qddot) + self.C.dot(self.qdot) + self.K.dot(self.q) - self.F
-
- return res
-
- def __TangentOperator(self):
- """ Description. """
-
- # The problem is linear, so the tangent operator is straightforward.
- St = self.betaPrime*self.M + self.gammaPrime*self.C + self.K
-
- return St
-
- def exit(self):
- """ Description. """
-
- print("\n**************** Exiting the structural tester solver ****************")
-
- def run(self,t0,t1):
- """ Description. """
-
- self.__temporalIteration()
-
- print("Time\tDisp 1\tDisp2\tVel 1\tVel2\tAcc 1\tAcc 2")
- print(str(t1) + '\t' + str(float(self.q[0])) + '\t' + str(float(self.q[1])) + '\t' + str(float(self.qdot[0])) + '\t' + str(float(self.qdot[1])) + '\t' + str(float(self.qddot[0])) + '\t' + str(float(self.qddot[1])))
-
- self.__computeInterfacePosVel(False)
-
- def setInitialDisplacements(self):
- """ Description. """
-
- self.__computeInterfacePosVel(True)
-
- def writeSolution(self, time, FSIIter, TimeIter, NbTimeIter):
- """ Description. """
-
- if time == 0:
- histFile = open('StructHistory.dat', "w")
- histFile.write("Time\tDisp 1\tDisp2\tVel 1\tVel2\tAcc 1\tAcc 2\tAccVar 1\tAccVar 2\n")
- else:
- histFile = open('StructHistory.dat', "a")
- histFile.write(str(time) + '\t' + str(float(self.q[0])) + '\t' + str(float(self.q[1])) + '\t' + str(float(self.qdot[0])) + '\t' + str(float(self.qdot[1])) + '\t' + str(float(self.qddot[0])) + '\t' + str(float(self.qddot[1])) + '\t' + str(float(self.a[0])) + '\t' + str(float(self.a[1])) + '\n')
- histFile.close()
-
- def updateSolution(self):
- """ Description. """
-
- self.q_n = np.copy(self.q)
- self.qdot_n = np.copy(self.qdot)
- self.qddot_n = np.copy(self.qddot)
- self.a_n = np.copy(self.a)
- self.__reset(self.q)
- self.__reset(self.qdot)
- self.__reset(self.qddot)
- self.__reset(self.a)
-
- makerID = self.markers.keys()[0]
- nodeList = self.markers[makerID]
-
- for iPoint in nodeList:
- self.node[iPoint].updateCoordVel()
-
- self.centerOfRotation_n = np.copy(self.centerOfRotation)
-
- def applyload(self, iVertex, fx, fy, fz, time):
- """ Description """
-
- makerID = self.markers.keys()[0]
- iPoint = self.getInterfaceNodeGlobalIndex(makerID, iVertex)
- self.node[iPoint].SetForce((fx,fy,fz))
-
- def getFSIMarkerID(self):
- """ Description. """
-
- list = self.markers.keys()
- return list[0]
-
- def getNumberOfSolidInterfaceNodes(self, markerID):
- """ Description. """
-
- return len(self.markers[markerID])
-
- def getInterfaceNodeGlobalIndex(self, markerID, iVertex):
- """ Description. """
-
- return self.markers[markerID][iVertex]
-
- def getInterfaceNodePosX(self, markerID, iVertex):
- """ Desciption. """
-
- iPoint = self.markers[markerID][iVertex]
- Coord = self.node[iPoint].GetCoord()
- return float(Coord[0])
-
- def getInterfaceNodePosY(self, markerID, iVertex):
- """ Desciption. """
-
- iPoint = self.markers[markerID][iVertex]
- Coord = self.node[iPoint].GetCoord()
- return float(Coord[1])
-
- def getInterfaceNodePosZ(self, markerID, iVertex):
- """ Desciption. """
-
- iPoint = self.markers[markerID][iVertex]
- Coord = self.node[iPoint].GetCoord()
- return float(Coord[2])
-
- def getInterfaceNodeDispX(self, markerID, iVertex):
- """ Desciption. """
-
- iPoint = self.markers[markerID][iVertex]
- Coord = self.node[iPoint].GetCoord()
- Coord0 = self.node[iPoint].GetCoord0()
- return float(Coord[0]-Coord0[0])
-
- def getInterfaceNodeDispY(self, markerID, iVertex):
- """ Desciption. """
-
- iPoint = self.markers[markerID][iVertex]
- Coord = self.node[iPoint].GetCoord()
- Coord0 = self.node[iPoint].GetCoord0()
- return float(Coord[1]-Coord0[1])
-
- def getInterfaceNodeDispZ(self, markerID, iVertex):
- """ Desciption. """
-
- iPoint = self.markers[markerID][iVertex]
- Coord = self.node[iPoint].GetCoord()
- Coord0 = self.node[iPoint].GetCoord0()
- return float(Coord[2]-Coord0[2])
-
- def getInterfaceNodeVelX(self, markerID, iVertex):
- """ Description """
-
- iPoint = self.markers[markerID][iVertex]
- Vel = self.node[iPoint].GetVel()
- return float(Vel[0])
-
- def getInterfaceNodeVelY(self, markerID, iVertex):
- """ Description """
-
- iPoint = self.markers[markerID][iVertex]
- Vel = self.node[iPoint].GetVel()
- return float(Vel[1])
-
- def getInterfaceNodeVelZ(self, markerID, iVertex):
- """ Description """
-
- iPoint = self.markers[markerID][iVertex]
- Vel = self.node[iPoint].GetVel()
- return float(Vel[2])
-
- def getInterfaceNodeVelXNm1(self, markerID, iVertex):
- """ Description """
-
- iPoint = self.markers[markerID][iVertex]
- Vel = self.node[iPoint].GetVel_n()
- return float(Vel[0])
-
- def getInterfaceNodeVelYNm1(self, markerID, iVertex):
- """ Description """
-
- iPoint = self.markers[markerID][iVertex]
- Vel = self.node[iPoint].GetVel_n()
- return float(Vel[1])
-
- def getInterfaceNodeVelZNm1(self, markerID, iVertex):
- """ Description """
-
- iPoint = self.markers[markerID][iVertex]
- Vel = self.node[iPoint].GetVel_n()
- return float(Vel[2])
-
- def getRotationCenterPosX(self):
- """ Description. """
-
- return float(self.centerOfRotation[0])
-
- def getRotationCenterPosY(self):
- """ Description. """
-
- return float(self.centerOfRotation[1])
-
- def getRotationCenterPosZ(self):
- """ Description. """
-
- return float(self.centerOfRotation[2])
diff --git a/SU2_PY/FSI/__init__.py b/SU2_PY/FSI/__init__.py
deleted file mode 100644
index f6e5483d7fcb..000000000000
--- a/SU2_PY/FSI/__init__.py
+++ /dev/null
@@ -1,4 +0,0 @@
-import io
-import util
-from FSIInterface import Interface
-import PitchPlungeAirfoilStructuralTester
diff --git a/SU2_PY/FSI/io/__init__.py b/SU2_PY/FSI/io/__init__.py
deleted file mode 100644
index b71c1c2d3ac0..000000000000
--- a/SU2_PY/FSI/io/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-from FSI_config import FSIConfig
diff --git a/SU2_PY/FSI/util/__init__.py b/SU2_PY/FSI/util/__init__.py
deleted file mode 100644
index 07a191da6d7b..000000000000
--- a/SU2_PY/FSI/util/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-from switch import switch
diff --git a/SU2_PY/FSI/FSIInterface.py b/SU2_PY/FSI_tools/FSIInterface.py
similarity index 65%
rename from SU2_PY/FSI/FSIInterface.py
rename to SU2_PY/FSI_tools/FSIInterface.py
index 3adaecc9df6d..951460377d6b 100644
--- a/SU2_PY/FSI/FSIInterface.py
+++ b/SU2_PY/FSI_tools/FSIInterface.py
@@ -2,12 +2,12 @@
## \file FSIInterface.py
# \brief FSI interface class that handles fluid/solid solvers synchronisation and communication.
-# \author David Thomas
+# \authors Nicola Fonzi, Vittorio Cavalieri based on the work of David Thomas
# \version 7.0.8 "Blackbird"
#
# SU2 Project Website: https://su2code.github.io
-#
-# The SU2 Project is maintained by the SU2 Foundation
+#
+# The SU2 Project is maintained by the SU2 Foundation
# (http://su2foundation.org)
#
# Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
@@ -16,7 +16,7 @@
# 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
@@ -42,16 +42,16 @@
# ----------------------------------------------------------------------
class Interface:
- """
+ """
FSI interface class that handles fluid/solid solvers synchronisation and communication
"""
-
+
def __init__(self, FSI_config, FluidSolver, SolidSolver, have_MPI):
- """
- Class constructor. Declare some variables and do some screen outputs.
- """
-
- if have_MPI == True:
+ """
+ Class constructor. Declare some variables and do some screen outputs.
+ """
+
+ if have_MPI:
from mpi4py import MPI
self.MPI = MPI
self.comm = MPI.COMM_WORLD #MPI World communicator
@@ -62,9 +62,9 @@ def __init__(self, FSI_config, FluidSolver, SolidSolver, have_MPI):
self.have_MPI = False
myid = 0
- self.rootProcess = 0 #the root process is chosen to be MPI rank = 0
+ self.rootProcess = 0 #the root process is chosen to be MPI rank = 0
- self.nDim = FSI_config['NDIM'] #problem dimension
+ self.nDim = FSI_config['NDIM'] #problem dimension
self.haveFluidSolver = False #True if the fluid solver is initialized on the current rank
self.haveSolidSolver = False #True if the solid solver is initialized on the current rank
@@ -87,16 +87,16 @@ def __init__(self, FSI_config, FluidSolver, SolidSolver, have_MPI):
self.SolidHaloNodeList = {} #contains the the indices (solid solver indexing) of the halo nodes for each partition
self.solidIndexing = {} #links between the solid solver indexing and the FSI indexing for the interface nodes
- self.nLocalFluidInterfaceNodes = 0 #number of nodes (halo nodes included) on the fluid interface, on each partition
+ self.nLocalFluidInterfaceNodes = 0 #number of nodes (halo nodes included) on the fluid interface, on each partition
self.nLocalFluidInterfaceHaloNode = 0 #number of halo nodes on the fluid intrface, on each partition
self.nLocalFluidInterfacePhysicalNodes = 0 #number of physical (= non halo) nodes on the fluid interface, on each partition
- self.nFluidInterfaceNodes = 0 #number of nodes on the fluid interface, sum over all the partitions
+ self.nFluidInterfaceNodes = 0 #number of nodes on the fluid interface, sum over all the partitions
self.nFluidInterfacePhysicalNodes = 0 #number of physical nodes on the fluid interface, sum over all partitions
self.nLocalSolidInterfaceNodes = 0 #number of physical nodes on the solid interface, on each partition
self.nLocalSolidInterfaceHaloNode = 0 #number of halo nodes on the solid intrface, on each partition
self.nLocalSolidInterfacePhysicalNodes = 0 #number of physical (= non halo) nodes on the solid interface, on each partition
- self.nSolidInterfaceNodes = 0 #number of nodes on the solid interface, sum over all partitions
+ self.nSolidInterfaceNodes = 0 #number of nodes on the solid interface, sum over all partitions
self.nSolidInterfacePhysicalNodes = 0 #number of physical nodes on the solid interface, sum over all partitions
if FSI_config['MATCHING_MESH'] == 'NO' and (FSI_config['MESH_INTERP_METHOD'] == 'RBF' or FSI_config['MESH_INTERP_METHOD'] == 'TPS'):
@@ -127,7 +127,7 @@ def __init__(self, FSI_config, FluidSolver, SolidSolver, have_MPI):
self.solidInterfaceResidualnM1_array_X = None #solid interface position residual at the previous BGS iteration
self.solidInterfaceResidualnM1_array_Y = None
self.solidInterfaceResidualnM1_array_Z = None
-
+
self.fluidInterface_array_DispX = None #fluid interface displacement
self.fluidInterface_array_DispY = None
self.fluidInterface_array_DispZ = None
@@ -140,53 +140,57 @@ def __init__(self, FSI_config, FluidSolver, SolidSolver, have_MPI):
self.solidLoads_array_Y = None
self.solidLoads_array_Z = None
- self.aitkenParam = FSI_config['AITKEN_PARAM'] #relaxation parameter for the BGS method
- self.FSIIter = 0 #current FSI iteration
+ self.aitkenParam = FSI_config['AITKEN_PARAM'] #relaxation parameter for the BGS method
+ self.FSIIter = 0 #current FSI iteration
self.unsteady = False #flag for steady or unsteady simulation (default is steady)
+ if FSI_config['CSD_SOLVER']=='IMPOSED':
+ self.ImposedMotion = True
+ else:
+ self.ImposedMotion = False
- # ---Some screen output ---
- self.MPIPrint('Fluid solver : SU2_CFD')
- self.MPIPrint('Solid solver : {}'.format(FSI_config['CSD_SOLVER']))
+ # ---Some screen output ---
+ self.MPIPrint('Fluid solver : SU2_CFD')
+ self.MPIPrint('Solid solver : {}'.format(FSI_config['CSD_SOLVER']))
- if FSI_config['TIME_MARCHING'] == 'YES':
+ if FSI_config['TIME_MARCHING'] == 'YES':
self.MPIPrint('Unsteady coupled simulation with physical time step : {} s'.format(FSI_config['UNST_TIMESTEP']))
self.unsteady = True
- else:
- self.MPIPrint('Steady coupled simulation')
+ else:
+ self.MPIPrint('Steady coupled simulation')
- if FSI_config['MATCHING_MESH'] == 'YES':
- self.MPIPrint('Matching fluid-solid interface')
- else:
+ if FSI_config['MATCHING_MESH'] == 'YES':
+ self.MPIPrint('Matching fluid-solid interface')
+ else:
if FSI_config['MESH_INTERP_METHOD'] == 'TPS':
- self.MPIPrint('Non matching fluid-solid interface with Thin Plate Spline interpolation')
+ self.MPIPrint('Non matching fluid-solid interface with Thin Plate Spline interpolation')
elif FSI_config['MESH_INTERP_METHOD'] == 'RBF':
self.MPIPrint('Non matching fluid-solid interface with Radial Basis Function interpolation')
self.RBF_rad = FSI_config['RBF_RADIUS']
- self.MPIPrint('Radius value : {}'.format(self.RBF_rad))
+ self.MPIPrint('Radius value : {}'.format(self.RBF_rad))
else:
- self.MPIPrint('Non matching fluid-solid interface with Nearest Neighboor interpolation')
+ self.MPIPrint('Non matching fluid-solid interface with Nearest Neighboor interpolation')
- self.MPIPrint('Solid predictor : {}'.format(FSI_config['DISP_PRED']))
+ self.MPIPrint('Solid predictor : {}'.format(FSI_config['DISP_PRED']))
- self.MPIPrint('Maximum number of FSI iterations : {}'.format(FSI_config['NB_FSI_ITER']))
+ self.MPIPrint('Maximum number of FSI iterations : {}'.format(FSI_config['NB_FSI_ITER']))
- self.MPIPrint('FSI tolerance : {}'.format(FSI_config['FSI_TOLERANCE']))
+ self.MPIPrint('FSI tolerance : {}'.format(FSI_config['FSI_TOLERANCE']))
- if FSI_config['AITKEN_RELAX'] == 'STATIC':
- self.MPIPrint('Static Aitken under-relaxation with constant parameter {}'.format(FSI_config['AITKEN_PARAM']))
- elif FSI_config['AITKEN_RELAX'] == 'DYNAMIC':
- self.MPIPrint('Dynamic Aitken under-relaxation with initial parameter {}'.format(FSI_config['AITKEN_PARAM']))
- else:
- self.MPIPrint('No Aitken under-relaxation')
+ if FSI_config['AITKEN_RELAX'] == 'STATIC':
+ self.MPIPrint('Static Aitken under-relaxation with constant parameter {}'.format(FSI_config['AITKEN_PARAM']))
+ elif FSI_config['AITKEN_RELAX'] == 'DYNAMIC':
+ self.MPIPrint('Dynamic Aitken under-relaxation with initial parameter {}'.format(FSI_config['AITKEN_PARAM']))
+ else:
+ self.MPIPrint('No Aitken under-relaxation')
self.MPIPrint('FSI interface is set')
def MPIPrint(self, message):
- """
+ """
Print a message on screen only from the master process.
"""
- if self.have_MPI == True:
+ if self.have_MPI:
myid = self.comm.Get_rank()
else:
myid = 0
@@ -198,86 +202,86 @@ def MPIBarrier(self):
"""
Perform a synchronization barrier in case of parallel run with MPI.
"""
-
- if self.have_MPI == True:
+
+ if self.have_MPI:
self.comm.barrier()
def connect(self, FSI_config, FluidSolver, SolidSolver):
- """
- Connection between solvers.
- Creates the communication support between the two solvers.
- Gets information about f/s interfaces from the two solvers.
- """
- if self.have_MPI == True:
+ """
+ Connection between solvers.
+ Creates the communication support between the two solvers.
+ Gets information about f/s interfaces from the two solvers.
+ """
+ if self.have_MPI:
myid = self.comm.Get_rank()
- MPIsize = self.comm.Get_size()
+ MPIsize = self.comm.Get_size()
else:
myid = 0
MPIsize = 1
-
- # --- Identify the fluid and solid interfaces and store the number of nodes on both sides (and for each partition) ---
+
+ # --- Identify the fluid and solid interfaces and store the number of nodes on both sides (and for each partition) ---
self.fluidInterfaceIdentifier = None
self.nLocalFluidInterfaceNodes = 0
- if FluidSolver != None:
- print('Fluid solver is initialized on process {}'.format(myid))
+ if FluidSolver is not None:
+ print('Fluid solver is initialized on process {}'.format(myid))
self.haveFluidSolver = True
- allMovingMarkersTags = FluidSolver.GetAllMovingMarkersTag()
+ allMovingMarkersTags = FluidSolver.GetAllDeformMeshMarkersTag()
allMarkersID = FluidSolver.GetAllBoundaryMarkers()
if not allMovingMarkersTags:
raise Exception('No interface for FSI was defined.')
else:
if allMovingMarkersTags[0] in allMarkersID.keys():
self.fluidInterfaceIdentifier = allMarkersID[allMovingMarkersTags[0]]
- if self.fluidInterfaceIdentifier != None:
- self.nLocalFluidInterfaceNodes = FluidSolver.GetNumberVertices(self.fluidInterfaceIdentifier)
- if self.nLocalFluidInterfaceNodes != 0:
+ if self.fluidInterfaceIdentifier is not None:
+ self.nLocalFluidInterfaceNodes = FluidSolver.GetNumberVertices(self.fluidInterfaceIdentifier)
+ if self.nLocalFluidInterfaceNodes != 0:
self.haveFluidInterface = True
- print('Number of interface fluid nodes (halo nodes included) on proccess {} : {}'.format(myid,self.nLocalFluidInterfaceNodes))
- else:
- pass
+ print('Number of interface fluid nodes (halo nodes included) on proccess {} : {}'.format(myid,self.nLocalFluidInterfaceNodes))
+ else:
+ pass
- if SolidSolver != None:
- print('Solid solver is initialized on process {}'.format(myid))
+ if SolidSolver is not None:
+ print('Solid solver is initialized on process {}'.format(myid))
self.haveSolidSolver = True
- self.solidInterfaceIdentifier = SolidSolver.getFSIMarkerID()
- self.nLocalSolidInterfaceNodes = SolidSolver.getNumberOfSolidInterfaceNodes(self.solidInterfaceIdentifier)
- if self.nLocalSolidInterfaceNodes != 0:
+ self.solidInterfaceIdentifier = SolidSolver.getFSIMarkerID()
+ self.nLocalSolidInterfaceNodes = SolidSolver.getNumberOfSolidInterfaceNodes(self.solidInterfaceIdentifier)
+ if self.nLocalSolidInterfaceNodes != 0:
self.haveSolidInterface = True
print('Number of interface solid nodes (halo nodes included) on proccess {} : {}'.format(myid,self.nLocalSolidInterfaceNodes))
- else:
- pass
+ else:
+ pass
# --- Exchange information about processors on which the solvers are defined and where the interface nodes are lying ---
- if self.have_MPI == True:
- if self.haveFluidSolver == True:
+ if self.have_MPI:
+ if self.haveFluidSolver:
sendBufFluid = np.array(int(1))
else:
sendBufFluid = np.array(int(0))
- if self.haveSolidSolver == True:
+ if self.haveSolidSolver:
sendBufSolid = np.array(int(1))
else:
sendBufSolid = np.array(int(0))
- if self.haveFluidInterface == True:
+ if self.haveFluidInterface:
sendBufFluidInterface = np.array(int(1))
else:
sendBufFluidInterface = np.array(int(0))
- if self.haveSolidInterface == True:
+ if self.haveSolidInterface:
sendBufSolidInterface = np.array(int(1))
else:
sendBufSolidInterface = np.array(int(0))
rcvBufFluid = np.zeros(MPIsize, dtype = int)
- rcvBufSolid = np.zeros(MPIsize, dtype = int)
+ rcvBufSolid = np.zeros(MPIsize, dtype = int)
rcvBufFluidInterface = np.zeros(MPIsize, dtype = int)
- rcvBufSolidInterface = np.zeros(MPIsize, dtype = int)
+ rcvBufSolidInterface = np.zeros(MPIsize, dtype = int)
self.comm.Allgather(sendBufFluid, rcvBufFluid)
self.comm.Allgather(sendBufSolid, rcvBufSolid)
self.comm.Allgather(sendBufFluidInterface, rcvBufFluidInterface)
self.comm.Allgather(sendBufSolidInterface, rcvBufSolidInterface)
for iProc in range(MPIsize):
- if rcvBufFluid[iProc] == 1:
+ if rcvBufFluid[iProc] == 1:
self.fluidSolverProcessors.append(iProc)
if rcvBufSolid[iProc] == 1:
- self.solidSolverProcessors.append(iProc)
+ self.solidSolverProcessors.append(iProc)
if rcvBufFluidInterface[iProc] == 1:
self.fluidInterfaceProcessors.append(iProc)
if rcvBufSolidInterface[iProc] == 1:
@@ -285,36 +289,30 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
del sendBufFluid, sendBufSolid, rcvBufFluid, rcvBufSolid, sendBufFluidInterface, sendBufSolidInterface, rcvBufFluidInterface, rcvBufSolidInterface
else:
self.fluidSolverProcessors.append(0)
- self.solidSolverProcessors.append(0)
+ self.solidSolverProcessors.append(0)
self.fluidInterfaceProcessors.append(0)
self.solidInterfaceProcessors.append(0)
- self.MPIBarrier()
-
- # --- Calculate the total number of nodes at the fluid interface (sum over all the partitions) ---
+ self.MPIBarrier()
+ # --- Calculate the total number of nodes at the fluid interface (sum over all the partitions) ---
# Calculate the number of halo nodes on each partition
self.nLocalFluidInterfaceHaloNode = 0
- for iVertex in range(self.nLocalFluidInterfaceNodes):
- if FluidSolver.IsAHaloNode(self.fluidInterfaceIdentifier, iVertex) == True:
+ for iVertex in range(self.nLocalFluidInterfaceNodes):
+ if FluidSolver.IsAHaloNode(self.fluidInterfaceIdentifier, iVertex):
GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex)
- self.FluidHaloNodeList[GlobalIndex] = iVertex
+ self.FluidHaloNodeList[GlobalIndex] = iVertex
self.nLocalFluidInterfaceHaloNode += 1
# Calculate the number of physical (= not halo) nodes on each partition
self.nLocalFluidInterfacePhysicalNodes = self.nLocalFluidInterfaceNodes - self.nLocalFluidInterfaceHaloNode
- if self.have_MPI == True:
+ if self.have_MPI:
self.FluidHaloNodeList = self.comm.allgather(self.FluidHaloNodeList)
else:
self.FluidHaloNodeList = [{}]
# Same thing for the solid part
self.nLocalSolidInterfaceHaloNode = 0
- #for iVertex in range(self.nLocalSolidInterfaceNodes):
- #if SoliddSolver.IsAHaloNode(self.fluidInterfaceIdentifier, iVertex) == True:
- #GlobalIndex = SolidSolver.GetVertexGlobalIndex(self.solidInterfaceIdentifier, iVertex)
- #self.SolidHaloNodeList[GlobalIndex] = iVertex
- #self.nLocalSolidInterfaceHaloNode += 1
self.nLocalSolidInterfacePhysicalNodes = self.nLocalSolidInterfaceNodes - self.nLocalSolidInterfaceHaloNode
- if self.have_MPI == True:
+ if self.have_MPI:
self.SolidHaloNodeList = self.comm.allgather(self.SolidHaloNodeList)
else:
self.SolidHaloNodeList = [{}]
@@ -323,11 +321,11 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
# --- Calculate the total number of nodes (with and without halo) at the fluid interface (sum over all the partitions) and broadcast the number accross all processors ---
sendBuffHalo = np.array(int(self.nLocalFluidInterfaceNodes))
sendBuffPhysical = np.array(int(self.nLocalFluidInterfacePhysicalNodes))
- rcvBuffHalo = np.zeros(1, dtype=int)
+ rcvBuffHalo = np.zeros(1, dtype=int)
rcvBuffPhysical = np.zeros(1, dtype=int)
- if self.have_MPI == True:
+ if self.have_MPI:
self.comm.barrier()
- self.comm.Allreduce(sendBuffHalo,rcvBuffHalo,op=self.MPI.SUM)
+ self.comm.Allreduce(sendBuffHalo,rcvBuffHalo,op=self.MPI.SUM)
self.comm.Allreduce(sendBuffPhysical,rcvBuffPhysical,op=self.MPI.SUM)
self.nFluidInterfaceNodes = rcvBuffHalo[0]
self.nFluidInterfacePhysicalNodes = rcvBuffPhysical[0]
@@ -339,11 +337,11 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
# Same thing for the solid part
sendBuffHalo = np.array(int(self.nLocalSolidInterfaceNodes))
sendBuffPhysical = np.array(int(self.nLocalSolidInterfacePhysicalNodes))
- rcvBuffHalo = np.zeros(1, dtype=int)
+ rcvBuffHalo = np.zeros(1, dtype=int)
rcvBuffPhysical = np.zeros(1, dtype=int)
- if self.have_MPI == True:
- self.comm.barrier()
- self.comm.Allreduce(sendBuffHalo,rcvBuffHalo,op=self.MPI.SUM)
+ if self.have_MPI:
+ self.comm.barrier()
+ self.comm.Allreduce(sendBuffHalo,rcvBuffHalo,op=self.MPI.SUM)
self.comm.Allreduce(sendBuffPhysical,rcvBuffPhysical,op=self.MPI.SUM)
self.nSolidInterfaceNodes = rcvBuffHalo[0]
self.nSolidInterfacePhysicalNodes = rcvBuffPhysical[0]
@@ -354,7 +352,7 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
# --- Store the number of physical interface nodes on each processor and allgather the information ---
self.fluidPhysicalInterfaceNodesDistribution = np.zeros(MPIsize, dtype=int)
- if self.have_MPI == True:
+ if self.have_MPI:
sendBuffPhysical = np.array(int(self.nLocalFluidInterfacePhysicalNodes))
self.comm.Allgather(sendBuffPhysical,self.fluidPhysicalInterfaceNodesDistribution)
del sendBuffPhysical
@@ -363,7 +361,7 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
# Same thing for the solid part
self.solidPhysicalInterfaceNodesDistribution = np.zeros(MPIsize, dtype=int)
- if self.have_MPI == True:
+ if self.have_MPI:
sendBuffPhysical = np.array(int(self.nLocalSolidInterfaceNodes))
self.comm.Allgather(sendBuffPhysical,self.solidPhysicalInterfaceNodesDistribution)
del sendBuffPhysical
@@ -371,11 +369,11 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
self.solidPhysicalInterfaceNodesDistribution[0] = self.nSolidInterfacePhysicalNodes
# --- Calculate and store the global indexing of interface physical nodes on each processor and allgather the information ---
- if self.have_MPI == True:
+ if self.have_MPI:
if myid in self.fluidInterfaceProcessors:
globalIndexStart = 0
for iProc in range(myid):
- globalIndexStart += self.fluidPhysicalInterfaceNodesDistribution[iProc]
+ globalIndexStart += self.fluidPhysicalInterfaceNodesDistribution[iProc]
globalIndexStop = globalIndexStart + self.nLocalFluidInterfacePhysicalNodes-1
else:
globalIndexStart = 0
@@ -387,9 +385,9 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
temp[0] = [0,self.nLocalFluidInterfacePhysicalNodes-1]
self.fluidGlobalIndexRange = list()
self.fluidGlobalIndexRange.append(temp)
-
- # Same thing for the solid part
- if self.have_MPI == True:
+
+ # Same thing for the solid part
+ if self.have_MPI:
if myid in self.solidInterfaceProcessors:
globalIndexStart = 0
for iProc in range(myid):
@@ -404,17 +402,17 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
temp = {}
temp[0] = [0,self.nSolidInterfacePhysicalNodes-1]
self.solidGlobalIndexRange = list()
- self.solidGlobalIndexRange.append(temp)
+ self.solidGlobalIndexRange.append(temp)
- self.MPIPrint('Total number of fluid interface nodes (halo nodes included) : {}'.format(self.nFluidInterfaceNodes))
- self.MPIPrint('Total number of solid interface nodes (halo nodes included) : {}'.format(self.nSolidInterfaceNodes))
+ self.MPIPrint('Total number of fluid interface nodes (halo nodes included) : {}'.format(self.nFluidInterfaceNodes))
+ self.MPIPrint('Total number of solid interface nodes (halo nodes included) : {}'.format(self.nSolidInterfaceNodes))
self.MPIPrint('Total number of fluid interface nodes : {}'.format(self.nFluidInterfacePhysicalNodes))
self.MPIPrint('Total number of solid interface nodes : {}'.format(self.nSolidInterfacePhysicalNodes))
- self.MPIBarrier()
+ self.MPIBarrier()
# --- Create all the PETSc vectors required for parallel communication and parallel mesh mapping/interpolation (working for serial too) ---
- if self.have_MPI == True:
+ if self.have_MPI:
self.solidInterface_array_DispX = PETSc.Vec().create(self.comm)
self.solidInterface_array_DispY = PETSc.Vec().create(self.comm)
self.solidInterface_array_DispZ = PETSc.Vec().create(self.comm)
@@ -432,10 +430,10 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
self.solidInterface_array_DispY.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
self.solidInterface_array_DispZ.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
self.solidInterface_array_DispX.set(0.0)
- self.solidInterface_array_DispY.set(0.0)
- self.solidInterface_array_DispZ.set(0.0)
+ self.solidInterface_array_DispY.set(0.0)
+ self.solidInterface_array_DispZ.set(0.0)
- if self.have_MPI == True:
+ if self.have_MPI:
self.fluidInterface_array_DispX = PETSc.Vec().create(self.comm)
self.fluidInterface_array_DispY = PETSc.Vec().create(self.comm)
self.fluidInterface_array_DispZ = PETSc.Vec().create(self.comm)
@@ -456,7 +454,7 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
self.fluidInterface_array_DispY.set(0.0)
self.fluidInterface_array_DispZ.set(0.0)
- if self.have_MPI == True:
+ if self.have_MPI:
self.fluidLoads_array_X = PETSc.Vec().create(self.comm)
self.fluidLoads_array_Y = PETSc.Vec().create(self.comm)
self.fluidLoads_array_Z = PETSc.Vec().create(self.comm)
@@ -473,8 +471,11 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
self.fluidLoads_array_X.setSizes(self.nFluidInterfacePhysicalNodes)
self.fluidLoads_array_Y.setSizes(self.nFluidInterfacePhysicalNodes)
self.fluidLoads_array_Z.setSizes(self.nFluidInterfacePhysicalNodes)
+ self.fluidLoads_array_X.set(0.0)
+ self.fluidLoads_array_Y.set(0.0)
+ self.fluidLoads_array_Z.set(0.0)
- if self.have_MPI == True:
+ if self.have_MPI:
self.solidLoads_array_X = PETSc.Vec().create(self.comm)
self.solidLoads_array_Y = PETSc.Vec().create(self.comm)
self.solidLoads_array_Z = PETSc.Vec().create(self.comm)
@@ -496,7 +497,7 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
self.solidLoads_array_Z.set(0.0)
# --- Create the PETSc vectors required for parallel relaxed BGS algo (working for serial too) ---
- if self.have_MPI == True:
+ if self.have_MPI:
self.solidInterfaceResidual_array_X = PETSc.Vec().create(self.comm)
self.solidInterfaceResidual_array_Y = PETSc.Vec().create(self.comm)
self.solidInterfaceResidual_array_Z = PETSc.Vec().create(self.comm)
@@ -513,8 +514,11 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
self.solidInterfaceResidual_array_X.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
self.solidInterfaceResidual_array_Y.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
self.solidInterfaceResidual_array_Z.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
+ self.solidInterfaceResidual_array_X.set(0.0)
+ self.solidInterfaceResidual_array_Y.set(0.0)
+ self.solidInterfaceResidual_array_Z.set(0.0)
- if self.have_MPI == True:
+ if self.have_MPI:
self.solidInterfaceResidualnM1_array_X = PETSc.Vec().create(self.comm)
self.solidInterfaceResidualnM1_array_Y = PETSc.Vec().create(self.comm)
self.solidInterfaceResidualnM1_array_Z = PETSc.Vec().create(self.comm)
@@ -536,30 +540,33 @@ def connect(self, FSI_config, FluidSolver, SolidSolver):
self.solidInterfaceResidualnM1_array_Z.set(0.0)
def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
- """
- Creates the one-to-one mapping between interfaces in case of matching meshes.
- Creates the interpolation rules between interfaces in case of non-matching meshes.
- """
- if self.have_MPI == True:
+ """
+ Creates the one-to-one mapping between interfaces in case of matching meshes.
+ Creates the interpolation rules between interfaces in case of non-matching meshes.
+ """
+ if self.have_MPI:
myid = self.comm.Get_rank()
- MPIsize = self.comm.Get_size()
+ MPIsize = self.comm.Get_size()
else:
myid = 0
MPIsize = 1
- # --- Get the fluid interface from fluid solver on each partition ---
- GlobalIndex = int()
+ # --- Get the fluid interface from fluid solver on each partition ---
+ GlobalIndex = int()
localIndex = 0
fluidIndexing_temp = {}
self.localFluidInterface_array_X_init = np.zeros((self.nLocalFluidInterfacePhysicalNodes))
self.localFluidInterface_array_Y_init = np.zeros((self.nLocalFluidInterfacePhysicalNodes))
self.localFluidInterface_array_Z_init = np.zeros((self.nLocalFluidInterfacePhysicalNodes))
for iVertex in range(self.nLocalFluidInterfaceNodes):
- GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex)
- posx = FluidSolver.GetVertexCoordX(self.fluidInterfaceIdentifier, iVertex)
- posy = FluidSolver.GetVertexCoordY(self.fluidInterfaceIdentifier, iVertex)
- posz = FluidSolver.GetVertexCoordZ(self.fluidInterfaceIdentifier, iVertex)
- if GlobalIndex in self.FluidHaloNodeList[myid].keys():
+ # Note that the fluid solver is separated in more processors outside the python script
+ # thus when, from a core, we request for the vertices on the interface, we only obtain
+ # those in that node
+ GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex)
+ posx = FluidSolver.GetVertexCoordX(self.fluidInterfaceIdentifier, iVertex)
+ posy = FluidSolver.GetVertexCoordY(self.fluidInterfaceIdentifier, iVertex)
+ posz = FluidSolver.GetVertexCoordZ(self.fluidInterfaceIdentifier, iVertex)
+ if GlobalIndex in self.FluidHaloNodeList[myid].keys():
self.haloNodesPositionsInit[GlobalIndex] = (posx, posy, posz)
else:
fluidIndexing_temp[GlobalIndex] = self.__getGlobalIndex('fluid', myid, localIndex)
@@ -567,26 +574,25 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
self.localFluidInterface_array_Y_init[localIndex] = posy
self.localFluidInterface_array_Z_init[localIndex] = posz
localIndex += 1
- if self.have_MPI == True:
+ if self.have_MPI:
fluidIndexing_temp = self.comm.allgather(fluidIndexing_temp)
for ii in range(len(fluidIndexing_temp)):
for key, value in fluidIndexing_temp[ii].items():
+ # This contains the link between the global index in python and that in SU2
self.fluidIndexing[key] = value
else:
self.fluidIndexing = fluidIndexing_temp.copy()
del fluidIndexing_temp
- # --- Get the solid interface from solid solver on each partition ---
+ # --- Get the solid interface from solid solver on each partition ---
localIndex = 0
solidIndexing_temp = {}
- self.localSolidInterface_array_X = np.zeros(self.nLocalSolidInterfaceNodes)
+ self.localSolidInterface_array_X = np.zeros(self.nLocalSolidInterfaceNodes)
self.localSolidInterface_array_Y = np.zeros(self.nLocalSolidInterfaceNodes)
self.localSolidInterface_array_Z = np.zeros(self.nLocalSolidInterfaceNodes)
for iVertex in range(self.nLocalSolidInterfaceNodes):
GlobalIndex = SolidSolver.getInterfaceNodeGlobalIndex(self.solidInterfaceIdentifier, iVertex)
- posx = SolidSolver.getInterfaceNodePosX(self.solidInterfaceIdentifier, iVertex)
- posy = SolidSolver.getInterfaceNodePosY(self.solidInterfaceIdentifier, iVertex)
- posz = SolidSolver.getInterfaceNodePosZ(self.solidInterfaceIdentifier, iVertex)
+ posx, posy, posz = SolidSolver.getInterfaceNodePos(self.solidInterfaceIdentifier, iVertex)
if GlobalIndex in self.SolidHaloNodeList[myid].keys():
pass
else:
@@ -595,7 +601,7 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
self.localSolidInterface_array_Y[localIndex] = posy
self.localSolidInterface_array_Z[localIndex] = posz
localIndex += 1
- if self.have_MPI == True:
+ if self.have_MPI:
solidIndexing_temp = self.comm.allgather(solidIndexing_temp)
for ii in range(len(solidIndexing_temp)):
for key, value in solidIndexing_temp[ii].items():
@@ -605,14 +611,14 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
del solidIndexing_temp
- # --- Create the PETSc parallel interpolation matrix ---
+ # --- Create the PETSc parallel interpolation matrix ---
if FSI_config['MATCHING_MESH'] == 'NO' and (FSI_config['MESH_INTERP_METHOD'] == 'RBF' or FSI_config['MESH_INTERP_METHOD'] == 'TPS'):
- if self.have_MPI == True:
+ if self.have_MPI:
self.MappingMatrixA = PETSc.Mat().create(self.comm)
self.MappingMatrixB = PETSc.Mat().create(self.comm)
self.MappingMatrixA_T = PETSc.Mat().create(self.comm)
self.MappingMatrixB_T = PETSc.Mat().create(self.comm)
- if FSI_config['MESH_INTERP_METHOD'] == 'RBF' :
+ if FSI_config['MESH_INTERP_METHOD'] == 'RBF' :
self.MappingMatrixA.setType('mpiaij')
self.MappingMatrixB.setType('mpiaij')
self.MappingMatrixA_T.setType('mpiaij')
@@ -627,7 +633,7 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
self.MappingMatrixB = PETSc.Mat().create()
self.MappingMatrixA_T = PETSc.Mat().create()
self.MappingMatrixB_T = PETSc.Mat().create()
- if FSI_config['MESH_INTERP_METHOD'] == 'RBF' :
+ if FSI_config['MESH_INTERP_METHOD'] == 'RBF' :
self.MappingMatrixA.setType('aij')
self.MappingMatrixB.setType('aij')
self.MappingMatrixA_T.setType('aij')
@@ -637,20 +643,20 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
self.MappingMatrixB.setType('aij')
self.MappingMatrixA_T.setType('aij')
self.MappingMatrixB_T.setType('aij')
- self.MappingMatrixA.setSizes((self.nSolidInterfacePhysicalNodes+self.d_RBF, self.nSolidInterfacePhysicalNodes+self.d_RBF))
+ self.MappingMatrixA.setSizes((self.nSolidInterfacePhysicalNodes+self.d_RBF, self.nSolidInterfacePhysicalNodes+self.d_RBF))
self.MappingMatrixA.setUp()
self.MappingMatrixA.setOption(PETSc.Mat().Option.NEW_NONZERO_ALLOCATION_ERR, False)
- self.MappingMatrixB.setSizes((self.nFluidInterfacePhysicalNodes, self.nSolidInterfacePhysicalNodes+self.d_RBF))
+ self.MappingMatrixB.setSizes((self.nFluidInterfacePhysicalNodes, self.nSolidInterfacePhysicalNodes+self.d_RBF))
self.MappingMatrixB.setUp()
self.MappingMatrixB.setOption(PETSc.Mat().Option.NEW_NONZERO_ALLOCATION_ERR, False)
- self.MappingMatrixA_T.setSizes((self.nSolidInterfacePhysicalNodes+self.d_RBF, self.nSolidInterfacePhysicalNodes+self.d_RBF))
+ self.MappingMatrixA_T.setSizes((self.nSolidInterfacePhysicalNodes+self.d_RBF, self.nSolidInterfacePhysicalNodes+self.d_RBF))
self.MappingMatrixA_T.setUp()
self.MappingMatrixA_T.setOption(PETSc.Mat().Option.NEW_NONZERO_ALLOCATION_ERR, False)
- self.MappingMatrixB_T.setSizes((self.nSolidInterfacePhysicalNodes+self.d_RBF, self.nFluidInterfacePhysicalNodes))
+ self.MappingMatrixB_T.setSizes((self.nSolidInterfacePhysicalNodes+self.d_RBF, self.nFluidInterfacePhysicalNodes))
self.MappingMatrixB_T.setUp()
self.MappingMatrixB_T.setOption(PETSc.Mat().Option.NEW_NONZERO_ALLOCATION_ERR, False)
else:
- if self.have_MPI == True:
+ if self.have_MPI:
self.MappingMatrix = PETSc.Mat().create(self.comm)
self.MappingMatrix_T = PETSc.Mat().create(self.comm)
self.MappingMatrix.setType('mpiaij')
@@ -660,32 +666,42 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
self.MappingMatrix_T = PETSc.Mat().create()
self.MappingMatrix.setType('aij')
self.MappingMatrix_T.setType('aij')
- self.MappingMatrix.setSizes((self.nFluidInterfacePhysicalNodes, self.nSolidInterfacePhysicalNodes))
+ self.MappingMatrix.setSizes((self.nFluidInterfacePhysicalNodes, self.nSolidInterfacePhysicalNodes))
self.MappingMatrix.setUp()
self.MappingMatrix.setOption(PETSc.Mat().Option.NEW_NONZERO_ALLOCATION_ERR, False)
- self.MappingMatrix_T.setSizes((self.nSolidInterfacePhysicalNodes, self.nFluidInterfacePhysicalNodes))
+ self.MappingMatrix_T.setSizes((self.nSolidInterfacePhysicalNodes, self.nFluidInterfacePhysicalNodes))
self.MappingMatrix_T.setUp()
self.MappingMatrix_T.setOption(PETSc.Mat().Option.NEW_NONZERO_ALLOCATION_ERR, False)
-
-
+
+
# --- Fill the interpolation matrix in parallel (working in serial too) ---
if FSI_config['MATCHING_MESH'] == 'NO' and (FSI_config['MESH_INTERP_METHOD'] == 'RBF' or FSI_config['MESH_INTERP_METHOD'] == 'TPS'):
self.MPIPrint('Building interpolation matrices...')
- if self.have_MPI == True:
+ if self.have_MPI:
for iProc in self.solidInterfaceProcessors:
if myid == iProc:
- for jProc in self.solidInterfaceProcessors:
- self.comm.Send(self.localSolidInterface_array_X, dest=jProc, tag=1)
- self.comm.Send(self.localSolidInterface_array_Y, dest=jProc, tag=2)
- self.comm.Send(self.localSolidInterface_array_Z, dest=jProc, tag=3)
+ for jProc in self.solidInterfaceProcessors:
+ if jProc != iProc:
+ self.comm.ssend(self.localSolidInterface_array_X, dest=jProc, tag=1)
+ self.comm.ssend(self.localSolidInterface_array_Y, dest=jProc, tag=2)
+ self.comm.ssend(self.localSolidInterface_array_Z, dest=jProc, tag=3)
+ else:
+ sizeOfBuff = self.solidPhysicalInterfaceNodesDistribution[iProc]
+ solidInterfaceBuffRcv_X = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_Y = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_Z = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_X = self.localSolidInterface_array_X
+ solidInterfaceBuffRcv_Y = self.localSolidInterface_array_Y
+ solidInterfaceBuffRcv_Z = self.localSolidInterface_array_Z
if myid in self.solidInterfaceProcessors:
- sizeOfBuff = self.solidPhysicalInterfaceNodesDistribution[iProc]
- solidInterfaceBuffRcv_X = np.zeros(sizeOfBuff)
- solidInterfaceBuffRcv_Y = np.zeros(sizeOfBuff)
- solidInterfaceBuffRcv_Z = np.zeros(sizeOfBuff)
- self.comm.Recv(solidInterfaceBuffRcv_X, iProc, tag=1)
- self.comm.Recv(solidInterfaceBuffRcv_Y, iProc, tag=2)
- self.comm.Recv(solidInterfaceBuffRcv_Z, iProc, tag=3)
+ if myid != iProc:
+ sizeOfBuff = self.solidPhysicalInterfaceNodesDistribution[iProc]
+ solidInterfaceBuffRcv_X = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_Y = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_Z = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_X = self.comm.recv(source=iProc, tag=1)
+ solidInterfaceBuffRcv_Y = self.comm.recv(source=iProc, tag=2)
+ solidInterfaceBuffRcv_Z = self.comm.recv(source=iProc, tag=3)
if FSI_config['MESH_INTERP_METHOD'] == 'RBF':
self.RBFMeshMapping_A(solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc, self.RBF_rad)
else:
@@ -703,22 +719,31 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
else:
self.MPIPrint("Building interpolation matrix...")
self.MPIBarrier()
-
- if self.have_MPI == True:
+ if self.have_MPI:
for iProc in self.solidInterfaceProcessors:
if myid == iProc:
for jProc in self.fluidInterfaceProcessors:
- self.comm.Send(self.localSolidInterface_array_X, dest=jProc, tag=1)
- self.comm.Send(self.localSolidInterface_array_Y, dest=jProc, tag=2)
- self.comm.Send(self.localSolidInterface_array_Z, dest=jProc, tag=3)
+ if jProc != iProc:
+ self.comm.ssend(self.localSolidInterface_array_X, dest=jProc, tag=1)
+ self.comm.ssend(self.localSolidInterface_array_Y, dest=jProc, tag=2)
+ self.comm.ssend(self.localSolidInterface_array_Z, dest=jProc, tag=3)
+ else:
+ sizeOfBuff = self.solidPhysicalInterfaceNodesDistribution[iProc]
+ solidInterfaceBuffRcv_X = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_Y = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_Z = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_X = self.localSolidInterface_array_X
+ solidInterfaceBuffRcv_Y = self.localSolidInterface_array_Y
+ solidInterfaceBuffRcv_Z = self.localSolidInterface_array_Z
if myid in self.fluidInterfaceProcessors:
- sizeOfBuff = self.solidPhysicalInterfaceNodesDistribution[iProc]
- solidInterfaceBuffRcv_X = np.zeros(sizeOfBuff)
- solidInterfaceBuffRcv_Y = np.zeros(sizeOfBuff)
- solidInterfaceBuffRcv_Z = np.zeros(sizeOfBuff)
- self.comm.Recv(solidInterfaceBuffRcv_X, iProc, tag=1)
- self.comm.Recv(solidInterfaceBuffRcv_Y, iProc, tag=2)
- self.comm.Recv(solidInterfaceBuffRcv_Z, iProc, tag=3)
+ if myid != iProc:
+ sizeOfBuff = self.solidPhysicalInterfaceNodesDistribution[iProc]
+ solidInterfaceBuffRcv_X = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_Y = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_Z = np.zeros(sizeOfBuff)
+ solidInterfaceBuffRcv_X = self.comm.recv(source=iProc, tag=1)
+ solidInterfaceBuffRcv_Y = self.comm.recv(source=iProc, tag=2)
+ solidInterfaceBuffRcv_Z = self.comm.recv(source=iProc, tag=3)
if FSI_config['MATCHING_MESH'] == 'NO':
if FSI_config['MESH_INTERP_METHOD'] == 'RBF':
self.RBFMeshMapping_B(solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc, self.RBF_rad)
@@ -726,7 +751,7 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
self.TPSMeshMapping_B(solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc)
else:
self.NearestNeighboorMeshMapping(solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc)
- else:
+ else:
self.matchingMeshMapping(solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc)
else:
if FSI_config['MATCHING_MESH'] == 'NO':
@@ -735,10 +760,10 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
elif FSI_config['MESH_INTERP_METHOD'] == 'TPS' :
self.TPSMeshMapping_B(self.localSolidInterface_array_X, self.localSolidInterface_array_Y, self.localSolidInterface_array_Z, 0)
else:
- self.NearestNeighboorMeshMapping(self.localSolidInterface_array_X, self.localSolidInterface_array_Y, self.localSolidInterface_array_Z, 0)
- else:
+ self.NearestNeighboorMeshMapping(self.localSolidInterface_array_X, self.localSolidInterface_array_Y, self.localSolidInterface_array_Z, 0)
+ else:
self.matchingMeshMapping(self.localSolidInterface_array_X, self.localSolidInterface_array_Y, self.localSolidInterface_array_Z, 0)
-
+
if FSI_config['MATCHING_MESH'] == 'NO' and (FSI_config['MESH_INTERP_METHOD'] == 'RBF' or FSI_config['MESH_INTERP_METHOD'] == 'TPS'):
self.MappingMatrixB.assemblyBegin()
self.MappingMatrixB.assemblyEnd()
@@ -751,37 +776,40 @@ def interfaceMapping(self,FluidSolver, SolidSolver, FSI_config):
self.MappingMatrix_T.assemblyBegin()
self.MappingMatrix_T.assemblyEnd()
self.MPIPrint("Interpolation matrix is built.")
-
+
self.MPIBarrier()
-
+
del self.localSolidInterface_array_X
del self.localSolidInterface_array_Y
del self.localSolidInterface_array_Z
+ del self.localFluidInterface_array_X_init
+ del self.localFluidInterface_array_Y_init
+ del self.localFluidInterface_array_Z_init
def matchingMeshMapping(self,solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc):
"""
Fill the mapping matrix in case of matching meshes at the f/s interface.
"""
- if self.have_MPI == True:
+ if self.have_MPI:
myid = self.comm.Get_rank()
else:
myid = 0
# --- Instantiate the spatial indexing ---
- prop_index = index.Property()
- prop_index.dimension = self.nDim
- SolidSpatialTree = index.Index(properties=prop_index)
-
+ prop_index = index.Property()
+ prop_index.dimension = self.nDim
+ SolidSpatialTree = index.Index(properties=prop_index)
+
nSolidNodes = solidInterfaceBuffRcv_X.shape[0]
for jVertex in range(nSolidNodes):
posX = solidInterfaceBuffRcv_X[jVertex]
posY = solidInterfaceBuffRcv_Y[jVertex]
posZ = solidInterfaceBuffRcv_Z[jVertex]
- if self.nDim == 2 :
- SolidSpatialTree.add(jVertex, (posX, posY))
- else :
- SolidSpatialTree.add(jVertex, (posX, posY, posZ))
+ if self.nDim == 2 :
+ SolidSpatialTree.add(jVertex, (posX, posY))
+ else :
+ SolidSpatialTree.add(jVertex, (posX, posY, posZ))
if self.nFluidInterfacePhysicalNodes != self.nSolidInterfacePhysicalNodes:
raise Exception("Fluid and solid interface must have the same number of nodes for matching meshes ! ")
@@ -807,35 +835,33 @@ def matchingMeshMapping(self,solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, s
self.MappingMatrix.setValue(iGlobalVertexFluid,jGlobalVertexSolid,1.0)
self.MappingMatrix_T.setValue(jGlobalVertexSolid, iGlobalVertexFluid,1.0)
- del solidInterfaceBuffRcv_X
- del solidInterfaceBuffRcv_Y
- del solidInterfaceBuffRcv_Z
-
def NearestNeighboorMeshMapping(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc):
"""
- Description
+ Interpolation based on the nearest neighboor.
+ For each node, the mesh is scanned to find the closed node to the first
+ one.
"""
- if self.have_MPI == True:
+ if self.have_MPI:
myid = self.comm.Get_rank()
else:
myid = 0
# --- Instantiate the spatial indexing ---
- prop_index = index.Property()
- prop_index.dimension = self.nDim
- SolidSpatialTree = index.Index(properties=prop_index)
-
+ prop_index = index.Property()
+ prop_index.dimension = self.nDim
+ SolidSpatialTree = index.Index(properties=prop_index)
+
nSolidNodes = solidInterfaceBuffRcv_X.shape[0]
for jVertex in range(nSolidNodes):
posX = solidInterfaceBuffRcv_X[jVertex]
posY = solidInterfaceBuffRcv_Y[jVertex]
posZ = solidInterfaceBuffRcv_Z[jVertex]
- if self.nDim == 2 :
- SolidSpatialTree.add(jVertex, (posX, posY))
- else :
- SolidSpatialTree.add(jVertex, (posX, posY, posZ))
+ if self.nDim == 2 :
+ SolidSpatialTree.add(jVertex, (posX, posY))
+ else :
+ SolidSpatialTree.add(jVertex, (posX, posY, posZ))
# --- For each fluid interface node, find the nearest solid interface node and fill the boolean mapping matrix ---
for iVertexFluid in range(self.nLocalFluidInterfacePhysicalNodes):
@@ -854,29 +880,31 @@ def NearestNeighboorMeshMapping(self, solidInterfaceBuffRcv_X, solidInterfaceBuf
def RBFMeshMapping_A(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc, rad):
"""
- Description
+ First part of the RBF mapping. This method provides the matrix required to
+ obtain, from the structural displacements, the loadings of the kernel
+ functions.
"""
- if self.have_MPI == True:
+ if self.have_MPI:
myid = self.comm.Get_rank()
else:
myid = 0
# --- Instantiate the spatial indexing ---
- prop_index = index.Property()
- prop_index.dimension = self.nDim
- SolidSpatialTree = index.Index(properties=prop_index)
-
+ prop_index = index.Property()
+ prop_index.dimension = self.nDim
+ SolidSpatialTree = index.Index(properties=prop_index)
+
nSolidNodes = solidInterfaceBuffRcv_X.shape[0]
for jVertex in range(nSolidNodes):
posX = solidInterfaceBuffRcv_X[jVertex]
posY = solidInterfaceBuffRcv_Y[jVertex]
posZ = solidInterfaceBuffRcv_Z[jVertex]
- if self.nDim == 2 :
- SolidSpatialTree.add(jVertex, (posX, posY))
- else :
- SolidSpatialTree.add(jVertex, (posX, posY, posZ))
+ if self.nDim == 2 :
+ SolidSpatialTree.add(jVertex, (posX, posY))
+ else :
+ SolidSpatialTree.add(jVertex, (posX, posY, posZ))
for iVertexSolid in range(self.nLocalSolidInterfaceNodes):
posX = self.localSolidInterface_array_X[iVertexSolid]
@@ -898,37 +926,40 @@ def RBFMeshMapping_A(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, sol
self.MappingMatrixA.setValue(iGlobalVertexSolid, nSolidNodes, 1.0)
self.MappingMatrixA.setValue(iGlobalVertexSolid, nSolidNodes+1, posX)
self.MappingMatrixA.setValue(iGlobalVertexSolid, nSolidNodes+2, posY)
- self.MappingMatrixA.setValue(iGlobalVertexSolid, nSolidNodes+3, posZ)
+ if self.nDim == 3:
+ self.MappingMatrixA.setValue(iGlobalVertexSolid, nSolidNodes+3, posZ)
self.MappingMatrixA_T.setValue(nSolidNodes, iGlobalVertexSolid, 1.0)
self.MappingMatrixA_T.setValue(nSolidNodes+1, iGlobalVertexSolid, posX)
self.MappingMatrixA_T.setValue(nSolidNodes+2, iGlobalVertexSolid, posY)
- self.MappingMatrixA_T.setValue(nSolidNodes+3, iGlobalVertexSolid, posZ)
+ if self.nDim == 3:
+ self.MappingMatrixA_T.setValue(nSolidNodes+3, iGlobalVertexSolid, posZ)
def RBFMeshMapping_B(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc, rad):
"""
- Description
+ Second part of the RBF mapping. This method provides the matrix required to
+ obtain, from the kernel function loadings, the fluid nodes displacements.
"""
- if self.have_MPI == True:
+ if self.have_MPI:
myid = self.comm.Get_rank()
else:
myid = 0
# --- Instantiate the spatial indexing ---
- prop_index = index.Property()
- prop_index.dimension = self.nDim
- SolidSpatialTree = index.Index(properties=prop_index)
-
+ prop_index = index.Property()
+ prop_index.dimension = self.nDim
+ SolidSpatialTree = index.Index(properties=prop_index)
+
nSolidNodes = solidInterfaceBuffRcv_X.shape[0]
for jVertex in range(nSolidNodes):
posX = solidInterfaceBuffRcv_X[jVertex]
posY = solidInterfaceBuffRcv_Y[jVertex]
posZ = solidInterfaceBuffRcv_Z[jVertex]
- if self.nDim == 2 :
- SolidSpatialTree.add(jVertex, (posX, posY))
- else :
- SolidSpatialTree.add(jVertex, (posX, posY, posZ))
+ if self.nDim == 2 :
+ SolidSpatialTree.add(jVertex, (posX, posY))
+ else :
+ SolidSpatialTree.add(jVertex, (posX, posY, posZ))
for iVertexFluid in range(self.nLocalFluidInterfacePhysicalNodes):
posX = self.localFluidInterface_array_X_init[iVertexFluid]
@@ -950,22 +981,26 @@ def RBFMeshMapping_B(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, sol
self.MappingMatrixB.setValue(iGlobalVertexFluid, nSolidNodes, 1.0)
self.MappingMatrixB.setValue(iGlobalVertexFluid, nSolidNodes+1, posX)
self.MappingMatrixB.setValue(iGlobalVertexFluid, nSolidNodes+2, posY)
- self.MappingMatrixB.setValue(iGlobalVertexFluid, nSolidNodes+3, posZ)
+ if self.nDim == 3:
+ self.MappingMatrixB.setValue(iGlobalVertexFluid, nSolidNodes+3, posZ)
self.MappingMatrixB_T.setValue(nSolidNodes, iGlobalVertexFluid, 1.0)
self.MappingMatrixB_T.setValue(nSolidNodes+1, iGlobalVertexFluid, posX)
self.MappingMatrixB_T.setValue(nSolidNodes+2, iGlobalVertexFluid, posY)
- self.MappingMatrixB_T.setValue(nSolidNodes+3, iGlobalVertexFluid, posZ)
+ if self.nDim == 3:
+ self.MappingMatrixB_T.setValue(nSolidNodes+3, iGlobalVertexFluid, posZ)
def TPSMeshMapping_A(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc):
"""
- Description
+ First part of the RBF mapping. This method provides the matrix required to
+ obtain, from the structural displacements, the loadings of the kernel
+ functions.
"""
- if self.have_MPI == True:
+ if self.have_MPI:
myid = self.comm.Get_rank()
else:
myid = 0
-
+
nSolidNodes = solidInterfaceBuffRcv_X.shape[0]
for iVertexSolid in range(self.nLocalSolidInterfaceNodes):
@@ -984,22 +1019,25 @@ def TPSMeshMapping_A(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, sol
self.MappingMatrixA.setValue(iGlobalVertexSolid, nSolidNodes, 1.0)
self.MappingMatrixA.setValue(iGlobalVertexSolid, nSolidNodes+1, posX)
self.MappingMatrixA.setValue(iGlobalVertexSolid, nSolidNodes+2, posY)
- self.MappingMatrixA.setValue(iGlobalVertexSolid, nSolidNodes+3, posZ)
+ if self.nDim == 3:
+ self.MappingMatrixA.setValue(iGlobalVertexSolid, nSolidNodes+3, posZ)
self.MappingMatrixA_T.setValue(nSolidNodes, iGlobalVertexSolid, 1.0)
self.MappingMatrixA_T.setValue(nSolidNodes+1, iGlobalVertexSolid, posX)
self.MappingMatrixA_T.setValue(nSolidNodes+2, iGlobalVertexSolid, posY)
- self.MappingMatrixA_T.setValue(nSolidNodes+3, iGlobalVertexSolid, posZ)
+ if self.nDim == 3:
+ self.MappingMatrixA_T.setValue(nSolidNodes+3, iGlobalVertexSolid, posZ)
def TPSMeshMapping_B(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, solidInterfaceBuffRcv_Z, iProc):
"""
- Description
+ Second part of the TPS mapping. This method provides the matrix required to
+ obtain, from the kernel function loadings, the fluid nodes displacements.
"""
- if self.have_MPI == True:
+ if self.have_MPI:
myid = self.comm.Get_rank()
else:
myid = 0
-
+
nSolidNodes = solidInterfaceBuffRcv_X.shape[0]
for iVertexFluid in range(self.nLocalFluidInterfacePhysicalNodes):
@@ -1018,20 +1056,23 @@ def TPSMeshMapping_B(self, solidInterfaceBuffRcv_X, solidInterfaceBuffRcv_Y, sol
self.MappingMatrixB.setValue(iGlobalVertexFluid, nSolidNodes, 1.0)
self.MappingMatrixB.setValue(iGlobalVertexFluid, nSolidNodes+1, posX)
self.MappingMatrixB.setValue(iGlobalVertexFluid, nSolidNodes+2, posY)
- self.MappingMatrixB.setValue(iGlobalVertexFluid, nSolidNodes+3, posZ)
+ if self.nDim == 3:
+ self.MappingMatrixB.setValue(iGlobalVertexFluid, nSolidNodes+3, posZ)
self.MappingMatrixB_T.setValue(nSolidNodes, iGlobalVertexFluid, 1.0)
self.MappingMatrixB_T.setValue(nSolidNodes+1, iGlobalVertexFluid, posX)
self.MappingMatrixB_T.setValue(nSolidNodes+2, iGlobalVertexFluid, posY)
- self.MappingMatrixB_T.setValue(nSolidNodes+3, iGlobalVertexFluid, posZ)
+ if self.nDim == 3:
+ self.MappingMatrixB_T.setValue(nSolidNodes+3, iGlobalVertexFluid, posZ)
def __CPC2(self, distance, rad):
"""
- Description.
+ This method provides the value of the kernel function given the euclidean
+ distance. The kernel function is the one used for RBF.
"""
phi = 0.0
eps = distance/rad
-
+
if eps < 1:
phi = ((1.0-eps)**4)*(4.0*eps+1.0)
else:
@@ -1041,23 +1082,24 @@ def __CPC2(self, distance, rad):
def __TPS(self, distance):
"""
- Description
+ This method provides the value of the kernel function given the euclidean
+ distance. The kernel function is the one used for TPS.
"""
phi = 0.0
-
+
if distance > 0.0:
phi = (distance**2)*np.log10(distance)
else:
phi = 0.0
- return phi
+ return phi
def interpolateSolidPositionOnFluidMesh(self, FSI_config):
- """
- Applies the one-to-one mapping or the interpolaiton rules from solid to fluid mesh.
- """
- if self.have_MPI == True:
+ """
+ Applies the one-to-one mapping or the interpolaiton rules from solid to fluid mesh.
+ """
+ if self.have_MPI:
myid = self.comm.Get_rank()
MPIsize = self.comm.Get_size()
else:
@@ -1067,7 +1109,7 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config):
# --- Interpolate (or map) in parallel the solid interface displacement on the fluid interface ---
if FSI_config['MATCHING_MESH'] == 'NO' and (FSI_config['MESH_INTERP_METHOD'] == 'RBF' or FSI_config['MESH_INTERP_METHOD'] == 'TPS'):
- if self.have_MPI == True:
+ if self.have_MPI:
gamma_array_DispX = PETSc.Vec().create(self.comm)
gamma_array_DispY = PETSc.Vec().create(self.comm)
gamma_array_DispZ = PETSc.Vec().create(self.comm)
@@ -1098,10 +1140,12 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config):
#print(KSP_solver.getInitialGuessNonzero())
KSP_solver.solve(self.solidInterface_array_DispX, gamma_array_DispX)
KSP_solver.solve(self.solidInterface_array_DispY, gamma_array_DispY)
- KSP_solver.solve(self.solidInterface_array_DispZ, gamma_array_DispZ)
+ if self.nDim==3:
+ KSP_solver.solve(self.solidInterface_array_DispZ, gamma_array_DispZ)
self.MappingMatrixB.mult(gamma_array_DispX, self.fluidInterface_array_DispX)
self.MappingMatrixB.mult(gamma_array_DispY, self.fluidInterface_array_DispY)
- self.MappingMatrixB.mult(gamma_array_DispZ, self.fluidInterface_array_DispZ)
+ if self.nDim==3:
+ self.MappingMatrixB.mult(gamma_array_DispZ, self.fluidInterface_array_DispZ)
gamma_array_DispX.destroy()
gamma_array_DispY.destroy()
gamma_array_DispZ.destroy()
@@ -1110,12 +1154,13 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config):
del gamma_array_DispY
del gamma_array_DispZ
del KSP_solver
- else:
+ else:
self.MappingMatrix.mult(self.solidInterface_array_DispX, self.fluidInterface_array_DispX)
self.MappingMatrix.mult(self.solidInterface_array_DispY, self.fluidInterface_array_DispY)
- self.MappingMatrix.mult(self.solidInterface_array_DispZ, self.fluidInterface_array_DispZ)
+ if self.nDim==3:
+ self.MappingMatrix.mult(self.solidInterface_array_DispZ, self.fluidInterface_array_DispZ)
- # --- Checking conservation ---
+ # --- Checking conservation ---
WSX = self.solidLoads_array_X.dot(self.solidInterface_array_DispX)
WSY = self.solidLoads_array_Y.dot(self.solidInterface_array_DispY)
WSZ = self.solidLoads_array_Z.dot(self.solidInterface_array_DispZ)
@@ -1124,14 +1169,16 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config):
WFY = self.fluidLoads_array_Y.dot(self.fluidInterface_array_DispY)
WFZ = self.fluidLoads_array_Z.dot(self.fluidInterface_array_DispZ)
- self.MPIPrint("Checking f/s interface conservation...")
- self.MPIPrint('Solid side (Wx, Wy, Wz) = ({}, {}, {})'.format(WSX, WSY, WSZ))
- self.MPIPrint('Fluid side (Wx, Wy, Wz) = ({}, {}, {})'.format(WFX, WFY, WFZ))
+ self.MPIPrint("Checking f/s interface conservation...")
+ self.MPIPrint('Solid side (Wx, Wy, Wz) = ({}, {}, {})'.format(WSX, WSY, WSZ))
+ self.MPIPrint('Fluid side (Wx, Wy, Wz) = ({}, {}, {})'.format(WFX, WFY, WFZ))
+
-
# --- Redistribute the interpolated fluid interface according to the partitions that own the fluid interface ---
# Gather the fluid interface on the master process
- if self.have_MPI == True:
+ # This is required because PETSc redistributes evenly in the cores, and does not use the same division
+ # of SU2, thus we need to redistribute
+ if self.have_MPI:
sendBuff_X = None
sendBuff_Y = None
sendBuff_Z = None
@@ -1156,7 +1203,7 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config):
displ = tuple(displ)
del sendBuffNumber, rcvBuffNumber
-
+
#print("DEBUG MESSAGE From proc {}, counts = {}".format(myid, counts))
#print("DEBUG MESSAGE From proc {}, displ = {}".format(myid, displ))
@@ -1176,19 +1223,26 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config):
sendBuff_Y[iVertex] = self.fluidInterface_array_DispY_recon[globalIndex]
sendBuff_Z[iVertex] = self.fluidInterface_array_DispZ_recon[globalIndex]
globalIndex += 1
- self.comm.Send(sendBuff_X, dest=iProc, tag = 1)
- self.comm.Send(sendBuff_Y, dest=iProc, tag = 2)
- self.comm.Send(sendBuff_Z, dest=iProc, tag = 3)
+ if iProc == self.rootProcess:
+ self.localFluidInterface_array_DispX = np.copy(sendBuff_X)
+ self.localFluidInterface_array_DispY = np.copy(sendBuff_Y)
+ self.localFluidInterface_array_DispZ = np.copy(sendBuff_Z)
+ if iProc != self.rootProcess:
+ self.comm.ssend(sendBuff_X, dest=iProc, tag = 1)
+ self.comm.ssend(sendBuff_Y, dest=iProc, tag = 2)
+ self.comm.ssend(sendBuff_Z, dest=iProc, tag = 3)
if myid in self.fluidInterfaceProcessors:
- self.localFluidInterface_array_DispX = np.zeros(self.nLocalFluidInterfacePhysicalNodes)
- self.localFluidInterface_array_DispY = np.zeros(self.nLocalFluidInterfacePhysicalNodes)
- self.localFluidInterface_array_DispZ = np.zeros(self.nLocalFluidInterfacePhysicalNodes)
- self.comm.Recv(self.localFluidInterface_array_DispX, source=self.rootProcess, tag = 1)
- self.comm.Recv(self.localFluidInterface_array_DispY, source=self.rootProcess, tag = 2)
- self.comm.Recv(self.localFluidInterface_array_DispZ, source=self.rootProcess, tag = 3)
+ if myid != self.rootProcess:
+ self.localFluidInterface_array_DispX = np.zeros(self.nLocalFluidInterfacePhysicalNodes)
+ self.localFluidInterface_array_DispY = np.zeros(self.nLocalFluidInterfacePhysicalNodes)
+ self.localFluidInterface_array_DispZ = np.zeros(self.nLocalFluidInterfacePhysicalNodes)
+ self.localFluidInterface_array_DispX = self.comm.recv(source=self.rootProcess, tag = 1)
+ self.localFluidInterface_array_DispY = self.comm.recv(source=self.rootProcess, tag = 2)
+ self.localFluidInterface_array_DispZ = self.comm.recv(source=self.rootProcess, tag = 3)
del sendBuff_X
del sendBuff_Y
del sendBuff_Z
+ self.comm.barrier()
else:
self.localFluidInterface_array_DispX = self.fluidInterface_array_DispX.getArray().copy()
self.localFluidInterface_array_DispY = self.fluidInterface_array_DispY.getArray().copy()
@@ -1197,7 +1251,7 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config):
# Special treatment for the halo nodes on the fluid interface
self.haloNodesDisplacements = {}
sendBuff = {}
- if self.have_MPI == True:
+ if self.have_MPI:
if myid == self.rootProcess:
for iProc in self.fluidInterfaceProcessors:
sendBuff = {}
@@ -1207,26 +1261,34 @@ def interpolateSolidPositionOnFluidMesh(self, FSI_config):
DispY = self.fluidInterface_array_DispY_recon[globalIndex]
DispZ = self.fluidInterface_array_DispZ_recon[globalIndex]
sendBuff[key] = (DispX, DispY, DispZ)
- self.comm.send(sendBuff, dest = iProc, tag=4)
+ if iProc == self.rootProcess:
+ self.haloNodesDisplacements = sendBuff
+ else:
+ self.comm.ssend(sendBuff, dest = iProc, tag=4)
if myid in self.fluidInterfaceProcessors:
- self.haloNodesDisplacements = self.comm.recv(source = self.rootProcess, tag = 4)
- del sendBuff
+ if myid != self.rootProcess:
+ self.haloNodesDisplacements = self.comm.recv(source = self.rootProcess, tag = 4)
+ self.comm.barrier()
+ del self.fluidInterface_array_DispX_recon
+ del self.fluidInterface_array_DispY_recon
+ del self.fluidInterface_array_DispZ_recon
+ del sendBuff
def interpolateFluidLoadsOnSolidMesh(self, FSI_config):
- """
- Applies the one-to-one mapping or the interpolaiton rules from fluid to solid mesh.
- """
- if self.have_MPI == True:
+ """
+ Applies the one-to-one mapping or the interpolaiton rules from fluid to solid mesh.
+ """
+ if self.have_MPI:
myid = self.comm.Get_rank()
MPIsize = self.comm.Get_size()
else:
myid = 0
MPIsize = 1
-
+
# --- Interpolate (or map) in parallel the fluid interface loads on the solid interface ---
- #self.MappingMatrix.transpose()
+ #self.MappingMatrix.transpose()
if FSI_config['MATCHING_MESH'] == 'NO' and (FSI_config['MESH_INTERP_METHOD'] == 'RBF' or FSI_config['MESH_INTERP_METHOD'] == 'TPS'):
- if self.have_MPI == True:
+ if self.have_MPI:
gamma_array_LoadX = PETSc.Vec().create(self.comm)
gamma_array_LoadY = PETSc.Vec().create(self.comm)
gamma_array_LoadZ = PETSc.Vec().create(self.comm)
@@ -1254,10 +1316,12 @@ def interpolateFluidLoadsOnSolidMesh(self, FSI_config):
KSP_solver.setFromOptions()
self.MappingMatrixB_T.mult(self.fluidLoads_array_X, gamma_array_LoadX)
self.MappingMatrixB_T.mult(self.fluidLoads_array_Y, gamma_array_LoadY)
- self.MappingMatrixB_T.mult(self.fluidLoads_array_Z, gamma_array_LoadZ)
+ if self.nDim==3:
+ self.MappingMatrixB_T.mult(self.fluidLoads_array_Z, gamma_array_LoadZ)
KSP_solver.solve(gamma_array_LoadX, self.solidLoads_array_X)
KSP_solver.solve(gamma_array_LoadY, self.solidLoads_array_Y)
- KSP_solver.solve(gamma_array_LoadZ, self.solidLoads_array_Z)
+ if self.nDim==3:
+ KSP_solver.solve(gamma_array_LoadZ, self.solidLoads_array_Z)
gamma_array_LoadX.destroy()
gamma_array_LoadY.destroy()
gamma_array_LoadZ.destroy()
@@ -1269,7 +1333,8 @@ def interpolateFluidLoadsOnSolidMesh(self, FSI_config):
else:
self.MappingMatrix_T.mult(self.fluidLoads_array_X, self.solidLoads_array_X)
self.MappingMatrix_T.mult(self.fluidLoads_array_Y, self.solidLoads_array_Y)
- self.MappingMatrix_T.mult(self.fluidLoads_array_Z, self.solidLoads_array_Z)
+ if self.nDim==3:
+ self.MappingMatrix_T.mult(self.fluidLoads_array_Z, self.solidLoads_array_Z)
# --- Redistribute the interpolated solid loads according to the partitions that own the solid interface ---
# Gather the solid loads on the master process
@@ -1280,10 +1345,10 @@ def interpolateFluidLoadsOnSolidMesh(self, FSI_config):
self.solidLoads_array_X_recon = None
self.solidLoads_array_Y_recon = None
self.solidLoads_array_Z_recon = None
- if myid == self.rootProcess:
- self.solidLoads_array_X_recon = np.zeros(self.nSolidInterfacePhysicalNodes+self.d_RBF)
- self.solidLoads_array_Y_recon = np.zeros(self.nSolidInterfacePhysicalNodes+self.d_RBF)
- self.solidLoads_array_Z_recon = np.zeros(self.nSolidInterfacePhysicalNodes+self.d_RBF)
+ if myid == self.rootProcess:
+ self.solidLoads_array_X_recon = np.zeros(self.nSolidInterfacePhysicalNodes+self.d_RBF)
+ self.solidLoads_array_Y_recon = np.zeros(self.nSolidInterfacePhysicalNodes+self.d_RBF)
+ self.solidLoads_array_Z_recon = np.zeros(self.nSolidInterfacePhysicalNodes+self.d_RBF)
myNumberOfNodes = self.solidLoads_array_X.getArray().shape[0]
sendBuffNumber = np.array([myNumberOfNodes], dtype=int)
rcvBuffNumber = np.zeros(MPIsize, dtype=int)
@@ -1293,9 +1358,9 @@ def interpolateFluidLoadsOnSolidMesh(self, FSI_config):
displ = np.zeros(MPIsize, dtype=int)
for ii in range(rcvBuffNumber.shape[0]):
displ[ii] = rcvBuffNumber[0:ii].sum()
- displ = tuple(displ)
+ displ = tuple(displ)
- del sendBuffNumber, rcvBuffNumber
+ del sendBuffNumber, rcvBuffNumber
self.comm.Gatherv(self.solidLoads_array_X.getArray(), [self.solidLoads_array_X_recon, counts, displ, self.MPI.DOUBLE], root=self.rootProcess)
self.comm.Gatherv(self.solidLoads_array_Y.getArray(), [self.solidLoads_array_Y_recon, counts, displ, self.MPI.DOUBLE], root=self.rootProcess)
@@ -1313,19 +1378,31 @@ def interpolateFluidLoadsOnSolidMesh(self, FSI_config):
sendBuff_Y[iVertex] = self.solidLoads_array_Y_recon[globalIndex]
sendBuff_Z[iVertex] = self.solidLoads_array_Z_recon[globalIndex]
globalIndex += 1
- self.comm.Send(sendBuff_X, dest=iProc, tag = 1)
- self.comm.Send(sendBuff_Y, dest=iProc, tag = 2)
- self.comm.Send(sendBuff_Z, dest=iProc, tag = 3)
+ if iProc != myid:
+ self.comm.ssend(sendBuff_X, dest=iProc, tag = 1)
+ self.comm.ssend(sendBuff_Y, dest=iProc, tag = 2)
+ self.comm.ssend(sendBuff_Z, dest=iProc, tag = 3)
+ else:
+ self.localSolidLoads_array_X = np.zeros(self.nLocalSolidInterfaceNodes)
+ self.localSolidLoads_array_Y = np.zeros(self.nLocalSolidInterfaceNodes)
+ self.localSolidLoads_array_Z = np.zeros(self.nLocalSolidInterfaceNodes)
+ self.localSolidLoads_array_X = sendBuff_X
+ self.localSolidLoads_array_Y = sendBuff_Y
+ self.localSolidLoads_array_Z = sendBuff_Z
if myid in self.solidInterfaceProcessors:
- self.localSolidLoads_array_X = np.zeros(self.nLocalSolidInterfaceNodes)
- self.localSolidLoads_array_Y = np.zeros(self.nLocalSolidInterfaceNodes)
- self.localSolidLoads_array_Z = np.zeros(self.nLocalSolidInterfaceNodes)
- self.comm.Recv(self.localSolidLoads_array_X, source=self.rootProcess, tag = 1)
- self.comm.Recv(self.localSolidLoads_array_Y, source=self.rootProcess, tag = 2)
- self.comm.Recv(self.localSolidLoads_array_Z, source=self.rootProcess, tag = 3)
+ if myid != self.rootProcess:
+ self.localSolidLoads_array_X = np.zeros(self.nLocalSolidInterfaceNodes)
+ self.localSolidLoads_array_Y = np.zeros(self.nLocalSolidInterfaceNodes)
+ self.localSolidLoads_array_Z = np.zeros(self.nLocalSolidInterfaceNodes)
+ self.localSolidLoads_array_X = self.comm.recv(source=self.rootProcess, tag = 1)
+ self.localSolidLoads_array_Y = self.comm.recv(source=self.rootProcess, tag = 2)
+ self.localSolidLoads_array_Z = self.comm.recv(source=self.rootProcess, tag = 3)
del sendBuff_X
del sendBuff_Y
del sendBuff_Z
+ del self.solidLoads_array_X_recon
+ del self.solidLoads_array_Y_recon
+ del self.solidLoads_array_Z_recon
else:
self.localSolidLoads_array_X = self.solidLoads_array_X.getArray().copy()
self.localSolidLoads_array_Y = self.solidLoads_array_Y.getArray().copy()
@@ -1334,66 +1411,24 @@ def interpolateFluidLoadsOnSolidMesh(self, FSI_config):
# Special treatment for the halo nodes on the fluid interface
# TODO when we will use parallel solid solver !!
-
- '''def getSolidInterfacePosition(self, SolidSolver):
- """
- Gets the current solid interface position from the solid solver.
- """
- if self.have_MPI == True:
- myid = self.comm.Get_rank()
- else:
- myid = 0
-
- # --- Get the solid interface position from the solid solver and directly fill the corresponding PETSc vector ---
- GlobalIndex = int()
- localIndex = 0
- for iVertex in range(self.nLocalSolidInterfaceNodes):
- GlobalIndex = SolidSolver.getInterfaceNodeGlobalIndex(self.solidInterfaceIdentifier, iVertex)
- if GlobalIndex in self.SolidHaloNodeList[myid].keys():
- pass
- else:
- newPosx = SolidSolver.getInterfaceNodePosX(self.solidInterfaceIdentifier, iVertex)
- newPosy = SolidSolver.getInterfaceNodePosY(self.solidInterfaceIdentifier, iVertex)
- newPosz = SolidSolver.getInterfaceNodePosZ(self.solidInterfaceIdentifier, iVertex)
- iGlobalVertex = self.__getGlobalIndex('solid', myid, localIndex)
- self.solidInterface_array_X.setValues([iGlobalVertex],newPosx)
- self.solidInterface_array_Y.setValues([iGlobalVertex],newPosy)
- self.solidInterface_array_Z.setValues([iGlobalVertex],newPosz)
- localIndex += 1
- #print("DEBUG MESSAGE From proc {} : In loop !".format(myid))
-
- #print("DEBUG MESSAGE From proc {} : Prepare for assembly !".format(myid))
-
- self.solidInterface_array_X.assemblyBegin()
- self.solidInterface_array_X.assemblyEnd()
- self.solidInterface_array_Y.assemblyBegin()
- self.solidInterface_array_Y.assemblyEnd()
- self.solidInterface_array_Z.assemblyBegin()
- self.solidInterface_array_Z.assemblyEnd()
-
- #print("DEBUG MESSAGE From PROC {} : Assembly is done !".format(myid))
- #print("DEBUG MESSAGE From PROC {} : array_X = {}".format(myid, self.solidInterface_array_X.getArray()))'''
-
def getSolidInterfaceDisplacement(self, SolidSolver):
- """
- Gets the current solid interface position from the solid solver.
- """
- if self.have_MPI == True:
- myid = self.comm.Get_rank()
+ """
+ Gets the current solid interface position from the solid solver.
+ """
+ if self.have_MPI:
+ myid = self.comm.Get_rank()
else:
myid = 0
-
+
# --- Get the solid interface position from the solid solver and directly fill the corresponding PETSc vector ---
GlobalIndex = int()
localIndex = 0
- for iVertex in range(self.nLocalSolidInterfaceNodes):
+ for iVertex in range(self.nLocalSolidInterfaceNodes):
GlobalIndex = SolidSolver.getInterfaceNodeGlobalIndex(self.solidInterfaceIdentifier, iVertex)
if GlobalIndex in self.SolidHaloNodeList[myid].keys():
pass
else:
- newDispx = SolidSolver.getInterfaceNodeDispX(self.solidInterfaceIdentifier, iVertex)
- newDispy = SolidSolver.getInterfaceNodeDispY(self.solidInterfaceIdentifier, iVertex)
- newDispz = SolidSolver.getInterfaceNodeDispZ(self.solidInterfaceIdentifier, iVertex)
+ newDispx, newDispy, newDispz = SolidSolver.getInterfaceNodeDisp(self.solidInterfaceIdentifier, iVertex)
iGlobalVertex = self.__getGlobalIndex('solid', myid, localIndex)
self.solidInterface_array_DispX.setValues([iGlobalVertex],newDispx)
self.solidInterface_array_DispY.setValues([iGlobalVertex],newDispy)
@@ -1408,10 +1443,10 @@ def getSolidInterfaceDisplacement(self, SolidSolver):
self.solidInterface_array_DispZ.assemblyEnd()
def getFluidInterfaceNodalForce(self, FSI_config, FluidSolver):
- """
- Gets the fluid interface loads from the fluid solver.
- """
- if self.have_MPI == True:
+ """
+ Gets the fluid interface loads from the fluid solver.
+ """
+ if self.have_MPI:
myid = self.comm.Get_rank()
else:
myid = 0
@@ -1422,27 +1457,20 @@ def getFluidInterfaceNodalForce(self, FSI_config, FluidSolver):
FZ = 0.0
# --- Get the fluid interface loads from the fluid solver and directly fill the corresponding PETSc vector ---
- for iVertex in range(self.nLocalFluidInterfaceNodes):
- halo = FluidSolver.ComputeVertexForces(self.fluidInterfaceIdentifier, iVertex) # !!we have to ignore halo node coming from mesh partitioning because they introduice non-physical forces
- if halo==False:
- if FSI_config['CSD_SOLVER'] == 'GETDP':
- newFx = FluidSolver.GetVertexForceDensityX(self.fluidInterfaceIdentifier, iVertex)
- newFy = FluidSolver.GetVertexForceDensityY(self.fluidInterfaceIdentifier, iVertex)
- newFz = FluidSolver.GetVertexForceDensityZ(self.fluidInterfaceIdentifier, iVertex)
- else:
- newFx = FluidSolver.GetVertexForceX(self.fluidInterfaceIdentifier, iVertex)
- newFy = FluidSolver.GetVertexForceY(self.fluidInterfaceIdentifier, iVertex)
- newFz = FluidSolver.GetVertexForceZ(self.fluidInterfaceIdentifier, iVertex)
- iGlobalVertex = self.__getGlobalIndex('fluid', myid, localIndex)
- self.fluidLoads_array_X.setValues([iGlobalVertex], newFx)
- self.fluidLoads_array_Y.setValues([iGlobalVertex], newFy)
- self.fluidLoads_array_Z.setValues([iGlobalVertex], newFz)
- FX += newFx
- FY += newFy
- FZ += newFz
- localIndex += 1
-
- if self.have_MPI == True:
+ for iVertex in range(self.nLocalFluidInterfaceNodes):
+ GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex)
+ if GlobalIndex not in self.FluidHaloNodeList[myid].keys():
+ loadX, loadY, loadZ = FluidSolver.GetFlowLoad(self.fluidInterfaceIdentifier, iVertex)
+ iGlobalVertex = self.__getGlobalIndex('fluid', myid, localIndex)
+ self.fluidLoads_array_X.setValues([iGlobalVertex], loadX)
+ self.fluidLoads_array_Y.setValues([iGlobalVertex], loadY)
+ self.fluidLoads_array_Z.setValues([iGlobalVertex], loadZ)
+ FX += loadX
+ FY += loadY
+ FZ += loadZ
+ localIndex += 1
+
+ if self.have_MPI:
FX = self.comm.allreduce(FX)
FY = self.comm.allreduce(FY)
FZ = self.comm.allreduce(FZ)
@@ -1457,79 +1485,68 @@ def getFluidInterfaceNodalForce(self, FSI_config, FluidSolver):
FX_b = self.fluidLoads_array_X.sum()
FY_b = self.fluidLoads_array_Y.sum()
FZ_b = self.fluidLoads_array_Z.sum()
-
+
def setFluidInterfaceVarCoord(self, FluidSolver):
- """
- Communicate the change of coordinates of the fluid interface to the fluid solver.
- Prepare the fluid solver for mesh deformation.
- """
- if self.have_MPI == True:
- myid = self.comm.Get_rank()
+ """
+ Communicate the change of coordinates of the fluid interface to the fluid solver.
+ Prepare the fluid solver for mesh deformation.
+ """
+ if self.have_MPI:
+ myid = self.comm.Get_rank()
else:
myid = 0
-
+
# --- Send the new fluid interface position to the fluid solver (on each partition, halo nodes included) ---
localIndex = 0
- for iVertex in range(self.nLocalFluidInterfaceNodes):
- GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex)
+ for iVertex in range(self.nLocalFluidInterfaceNodes):
+ GlobalIndex = FluidSolver.GetVertexGlobalIndex(self.fluidInterfaceIdentifier, iVertex)
if GlobalIndex in self.FluidHaloNodeList[myid].keys():
- posX0, posY0, posZ0 = self.haloNodesPositionsInit[GlobalIndex]
DispX, DispY, DispZ = self.haloNodesDisplacements[GlobalIndex]
- posX = posX0 + DispX
- posY = posY0 + DispY
- posZ = posZ0 + DispZ
- FluidSolver.SetVertexCoordX(self.fluidInterfaceIdentifier, iVertex, posX)
- FluidSolver.SetVertexCoordY(self.fluidInterfaceIdentifier, iVertex, posY)
- FluidSolver.SetVertexCoordZ(self.fluidInterfaceIdentifier, iVertex, posZ)
+ FluidSolver.SetMeshDisplacement(self.fluidInterfaceIdentifier, int(iVertex), DispX, DispY, DispZ)
else:
- posX = self.localFluidInterface_array_DispX[localIndex] + self.localFluidInterface_array_X_init[localIndex]
- posY = self.localFluidInterface_array_DispY[localIndex] + self.localFluidInterface_array_Y_init[localIndex]
- posZ = self.localFluidInterface_array_DispZ[localIndex] + self.localFluidInterface_array_Z_init[localIndex]
- FluidSolver.SetVertexCoordX(self.fluidInterfaceIdentifier, iVertex, posX)
- FluidSolver.SetVertexCoordY(self.fluidInterfaceIdentifier, iVertex, posY)
- FluidSolver.SetVertexCoordZ(self.fluidInterfaceIdentifier, iVertex, posZ)
+ DispX = self.localFluidInterface_array_DispX[localIndex]
+ DispY = self.localFluidInterface_array_DispY[localIndex]
+ DispZ = self.localFluidInterface_array_DispZ[localIndex]
+ FluidSolver.SetMeshDisplacement(self.fluidInterfaceIdentifier, int(iVertex), DispX, DispY, DispZ)
localIndex += 1
- # Prepares the mesh deformation in the fluid solver
- nodalVarCoordNorm = FluidSolver.SetVertexVarCoord(self.fluidInterfaceIdentifier, iVertex)
-
-
- def setSolidInterfaceLoads(self, SolidSolver, FSI_config, time):
- """
- Communicates the new solid interface loads to the solid solver.
- In case of rigid body motion, calculates the new resultant forces (lift, drag, ...).
- """
- if self.have_MPI == True:
- myid = self.comm.Get_rank()
+
+
+ def setSolidInterfaceLoads(self, SolidSolver, FSI_config):
+ """
+ Communicates the new solid interface loads to the solid solver.
+ In case of rigid body motion, calculates the new resultant forces (lift, drag, ...).
+ """
+ if self.have_MPI:
+ myid = self.comm.Get_rank()
else:
myid = 0
- FY = 0.0 # solid-side resultant forces
+ FY = 0.0 # solid-side resultant forces
FX = 0.0
FZ = 0.0
- FFX = 0.0 # fluid-side resultant forces
- FFY = 0.0
- FFZ = 0.0
+ FFX = 0.0 # fluid-side resultant forces
+ FFY = 0.0
+ FFZ = 0.0
# --- Check for total force conservation after interpolation
FFX = self.fluidLoads_array_X.sum()
FFY = self.fluidLoads_array_Y.sum()
FFZ = self.fluidLoads_array_Z.sum()
-
for iVertex in range(self.nLocalSolidInterfaceNodes):
FX += self.localSolidLoads_array_X[iVertex]
FY += self.localSolidLoads_array_Y[iVertex]
FZ += self.localSolidLoads_array_Z[iVertex]
- if self.have_MPI == True:
+ if self.have_MPI:
FX = self.comm.allreduce(FX)
FY = self.comm.allreduce(FY)
FZ = self.comm.allreduce(FZ)
- self.MPIPrint("Checking f/s interface total force...")
- self.MPIPrint('Solid side (Fx, Fy, Fz) = ({}, {}, {})'.format(FX, FY, FZ))
- self.MPIPrint('Fluid side (Fx, Fy, Fz) = ({}, {}, {})'.format(FFX, FFY, FFZ))
+ self.MPIPrint("Checking f/s interface total force...")
+ self.MPIPrint('Solid side (Fx, Fy, Fz) = ({}, {}, {})'.format(FX, FY, FZ))
+ self.MPIPrint('Fluid side (Fx, Fy, Fz) = ({}, {}, {})'.format(FFX, FFY, FFZ))
# --- Send the new solid interface loads to the solid solver (on each partition, halo nodes included) ---
GlobalIndex = int()
@@ -1541,28 +1558,25 @@ def setSolidInterfaceLoads(self, SolidSolver, FSI_config, time):
pass
else:
Fx = self.localSolidLoads_array_X[localIndex]
- Fy = self.localSolidLoads_array_Y[localIndex]
- Fz = self.localSolidLoads_array_Z[localIndex]
- SolidSolver.applyload(iVertex, Fx, Fy, Fz, time)
+ Fy = self.localSolidLoads_array_Y[localIndex]
+ Fz = self.localSolidLoads_array_Z[localIndex]
+ SolidSolver.applyload(iVertex, Fx, Fy, Fz)
localIndex += 1
- if FSI_config['CSD_SOLVER'] == 'NATIVE':
- SolidSolver.setGeneralisedForce()
- SolidSolver.setGeneralisedMoment()
def computeSolidInterfaceResidual(self, SolidSolver):
- """
- Computes the solid interface FSI displacement residual.
- """
+ """
+ Computes the solid interface FSI displacement residual.
+ """
- if self.have_MPI == True:
- myid = self.comm.Get_rank()
+ if self.have_MPI:
+ myid = self.comm.Get_rank()
else:
myid = 0
- normInterfaceResidualSquare = 0.0
+ normInterfaceResidualSquare = 0.0
# --- Create and fill the PETSc vector for the predicted solid interface position (predicted by the solid computation) ---
- if self.have_MPI == True:
+ if self.have_MPI:
predDisp_array_X = PETSc.Vec().create(self.comm)
predDisp_array_X.setType('mpi')
predDisp_array_Y = PETSc.Vec().create(self.comm)
@@ -1575,27 +1589,25 @@ def computeSolidInterfaceResidual(self, SolidSolver):
predDisp_array_Y = PETSc.Vec().create()
predDisp_array_Y.setType('seq')
predDisp_array_Z = PETSc.Vec().create()
- predDisp_array_Z.setType('seq')
+ predDisp_array_Z.setType('seq')
predDisp_array_X.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
predDisp_array_Y.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
predDisp_array_Z.setSizes(self.nSolidInterfacePhysicalNodes+self.d_RBF)
-
- if myid in self.solidSolverProcessors:
- for iVertex in range(self.nLocalSolidInterfaceNodes):
- predDispx = SolidSolver.getInterfaceNodeDispX(self.solidInterfaceIdentifier, iVertex)
- predDispy = SolidSolver.getInterfaceNodeDispY(self.solidInterfaceIdentifier, iVertex)
- predDispz = SolidSolver.getInterfaceNodeDispZ(self.solidInterfaceIdentifier, iVertex)
+
+ if myid in self.solidSolverProcessors:
+ for iVertex in range(self.nLocalSolidInterfaceNodes):
+ predDispx, predDispy, predDispz = SolidSolver.getInterfaceNodeDisp(self.solidInterfaceIdentifier, iVertex)
iGlobalVertex = self.__getGlobalIndex('solid', myid, iVertex)
predDisp_array_X.setValues([iGlobalVertex], predDispx)
predDisp_array_Y.setValues([iGlobalVertex], predDispy)
predDisp_array_Z.setValues([iGlobalVertex], predDispz)
-
- predDisp_array_X.assemblyBegin()
- predDisp_array_X.assemblyEnd()
- predDisp_array_Y.assemblyBegin()
- predDisp_array_Y.assemblyEnd()
- predDisp_array_Z.assemblyBegin()
- predDisp_array_Z.assemblyEnd()
+
+ predDisp_array_X.assemblyBegin()
+ predDisp_array_X.assemblyEnd()
+ predDisp_array_Y.assemblyBegin()
+ predDisp_array_Y.assemblyEnd()
+ predDisp_array_Z.assemblyBegin()
+ predDisp_array_Z.assemblyEnd()
# --- Calculate the residual (vector and norm) ---
self.solidInterfaceResidual_array_X = predDisp_array_X - self.solidInterface_array_DispX
@@ -1615,45 +1627,45 @@ def computeSolidInterfaceResidual(self, SolidSolver):
del predDisp_array_Y
del predDisp_array_Z
- return sqrt(normInterfaceResidualSquare)
+ return sqrt(normInterfaceResidualSquare)
def relaxSolidPosition(self,FSI_config):
- """
- Apply solid displacement under-relaxation.
- """
- if self.have_MPI == True:
- myid = self.comm.Get_rank()
+ """
+ Apply solid displacement under-relaxation.
+ """
+ if self.have_MPI:
+ myid = self.comm.Get_rank()
else:
myid = 0
# --- Set the Aitken coefficient for the relaxation ---
- if FSI_config['AITKEN_RELAX'] == 'STATIC':
- self.aitkenParam = FSI_config['AITKEN_PARAM']
- elif FSI_config['AITKEN_RELAX'] == 'DYNAMIC':
- self.setAitkenCoefficient(FSI_config)
- else:
- self.aitkenParam = 1.0
+ if FSI_config['AITKEN_RELAX'] == 'STATIC':
+ self.aitkenParam = FSI_config['AITKEN_PARAM']
+ elif FSI_config['AITKEN_RELAX'] == 'DYNAMIC':
+ self.setAitkenCoefficient(FSI_config)
+ else:
+ self.aitkenParam = 1.0
- self.MPIPrint('Aitken under-relaxation step with parameter {}'.format(self.aitkenParam))
+ self.MPIPrint('Aitken under-relaxation step with parameter {}'.format(self.aitkenParam))
# --- Relax the solid interface position ---
self.solidInterface_array_DispX += self.aitkenParam*self.solidInterfaceResidual_array_X
self.solidInterface_array_DispY += self.aitkenParam*self.solidInterfaceResidual_array_Y
self.solidInterface_array_DispZ += self.aitkenParam*self.solidInterfaceResidual_array_Z
-
+
def setAitkenCoefficient(self, FSI_config):
- """
- Computes the Aitken coefficients for solid displacement under-relaxation.
- """
+ """
+ Computes the Aitken coefficients for solid displacement under-relaxation.
+ """
+
+ deltaResNormSquare = 0.0
+ prodScalRes = 0.0
- deltaResNormSquare = 0.0
- prodScalRes = 0.0
-
# --- Create the PETSc vector for the difference between the residuals (current and previous FSI iter) ---
- if self.FSIIter == 0:
- self.aitkenParam = max(FSI_config['AITKEN_PARAM'], self.aitkenParam)
- else:
+ if self.FSIIter == 0:
+ self.aitkenParam = max(FSI_config['AITKEN_PARAM'], self.aitkenParam)
+ else:
if self.have_MPI:
deltaResx_array_X = PETSc.Vec().create(self.comm)
deltaResx_array_X.setType('mpi')
@@ -1688,9 +1700,9 @@ def setAitkenCoefficient(self, FSI_config):
deltaResNormSquare_X = (deltaResx_array_X.norm())**2
deltaResNormSquare_Y = (deltaResx_array_Y.norm())**2
deltaResNormSquare_Z = (deltaResx_array_Z.norm())**2
- deltaResNormSquare = deltaResNormSquare_X + deltaResNormSquare_Y + deltaResNormSquare_Z
+ deltaResNormSquare = deltaResNormSquare_X + deltaResNormSquare_Y + deltaResNormSquare_Z
- self.aitkenParam *= -prodScalRes/deltaResNormSquare
+ self.aitkenParam *= -prodScalRes/deltaResNormSquare
deltaResx_array_X.destroy()
deltaResx_array_Y.destroy()
@@ -1708,30 +1720,30 @@ def setAitkenCoefficient(self, FSI_config):
self.solidInterfaceResidual_array_Z.copy(self.solidInterfaceResidualnM1_array_Z)
def displacementPredictor(self, FSI_config , SolidSolver, deltaT):
- """
- Calculates a prediciton for the solid interface position for the next time step.
- """
+ """
+ Calculates a prediciton for the solid interface position for the next time step.
+ """
- if self.have_MPI == True:
- myid = self.comm.Get_rank()
+ if self.have_MPI:
+ myid = self.comm.Get_rank()
else:
myid = 0
- if FSI_config['DISP_PRED'] == 'FIRST_ORDER':
- self.MPIPrint("First order predictor")
- alpha_0 = 1.0
- alpha_1 = 0.0
- elif FSI_config['DISP_PRED'] == 'SECOND_ORDER':
- self.MPIPrint("Second order predictor")
- alpha_0 = 1.0
- alpha_1 = 0.5
- else:
- self.MPIPrint("No predictor")
- alpha_0 = 0.0
- alpha_1 = 0.0
+ if FSI_config['DISP_PRED'] == 'FIRST_ORDER':
+ self.MPIPrint("First order predictor")
+ alpha_0 = 1.0
+ alpha_1 = 0.0
+ elif FSI_config['DISP_PRED'] == 'SECOND_ORDER':
+ self.MPIPrint("Second order predictor")
+ alpha_0 = 1.0
+ alpha_1 = 0.5
+ else:
+ self.MPIPrint("No predictor")
+ alpha_0 = 0.0
+ alpha_1 = 0.0
# --- Create the PETSc vectors to store the solid interface velocity ---
- if self.have_MPI == True:
+ if self.have_MPI:
Vel_array_X = PETSc.Vec().create(self.comm)
Vel_array_X.setType('mpi')
Vel_array_Y = PETSc.Vec().create(self.comm)
@@ -1774,18 +1786,14 @@ def displacementPredictor(self, FSI_config , SolidSolver, deltaT):
# --- Fill the PETSc vectors ---
GlobalIndex = int()
localIndex = 0
- for iVertex in range(self.nLocalSolidInterfaceNodes):
- GlobalIndex = SolidSolver.getInterfaceNodeGlobalIndex(self.solidInterfaceIdentifier, iVertex)
+ for iVertex in range(self.nLocalSolidInterfaceNodes):
+ GlobalIndex = SolidSolver.getInterfaceNodeGlobalIndex(self.solidInterfaceIdentifier, iVertex)
if GlobalIndex in self.SolidHaloNodeList[myid].keys():
pass
else:
iGlobalVertex = self.__getGlobalIndex('solid', myid, localIndex)
- velx = SolidSolver.getInterfaceNodeVelX(self.solidInterfaceIdentifier, iVertex)
- vely = SolidSolver.getInterfaceNodeVelY(self.solidInterfaceIdentifier, iVertex)
- velz = SolidSolver.getInterfaceNodeVelZ(self.solidInterfaceIdentifier, iVertex)
- velxNm1 = SolidSolver.getInterfaceNodeVelXNm1(self.solidInterfaceIdentifier, iVertex)
- velyNm1 = SolidSolver.getInterfaceNodeVelYNm1(self.solidInterfaceIdentifier, iVertex)
- velzNm1 = SolidSolver.getInterfaceNodeVelZNm1(self.solidInterfaceIdentifier, iVertex)
+ velx, vely, velz = SolidSolver.getInterfaceNodeVel(self.solidInterfaceIdentifier, iVertex)
+ velxNm1, velyNm1, velzNm1 = SolidSolver.getInterfaceNodeVelNm1(self.solidInterfaceIdentifier, iVertex)
Vel_array_X.setValues([iGlobalVertex],velx)
Vel_array_Y.setValues([iGlobalVertex],vely)
Vel_array_Z.setValues([iGlobalVertex],velz)
@@ -1822,27 +1830,27 @@ def displacementPredictor(self, FSI_config , SolidSolver, deltaT):
del VelnM1_array_X, VelnM1_array_Y, VelnM1_array_Z
def writeFSIHistory(self, TimeIter, time, varCoordNorm, FSIConv):
- """
- Write the FSI history file of the computaion.
- """
+ """
+ Write the FSI history file of the computaion.
+ """
- if self.have_MPI == True:
+ if self.have_MPI:
myid = self.comm.Get_rank()
else:
myid = 0
-
+
if myid == self.rootProcess:
- if self.unsteady:
- if TimeIter == 0:
- histFile = open('FSIhistory.dat', "w")
+ if self.unsteady:
+ if TimeIter == 0:
+ histFile = open('FSIhistory.dat', "w")
histFile.write("TimeIter\tTime\tFSIRes\tFSINbIter\n")
- else:
- histFile = open('FSIhistory.dat', "a")
- if FSIConv:
- histFile.write(str(TimeIter) + '\t' + str(time) + '\t' + str(varCoordNorm) + '\t' + str(self.FSIIter+1) + '\n')
- else:
- histFile.write(str(TimeIter) + '\t' + str(time) + '\t' + str(varCoordNorm) + '\t' + str(self.FSIIter) + '\n')
- histFile.close()
+ else:
+ histFile = open('FSIhistory.dat', "a")
+ if FSIConv:
+ histFile.write(str(TimeIter) + '\t' + str(time) + '\t' + str(varCoordNorm) + '\t' + str(self.FSIIter+1) + '\n')
+ else:
+ histFile.write(str(TimeIter) + '\t' + str(time) + '\t' + str(varCoordNorm) + '\t' + str(self.FSIIter) + '\n')
+ histFile.close()
else:
if self.FSIIter == 0:
histFile = open('FSIhistory.dat', "w")
@@ -1851,13 +1859,16 @@ def writeFSIHistory(self, TimeIter, time, varCoordNorm, FSIConv):
histFile = open('FSIhistory.dat', "a")
histFile.write(str(self.FSIIter) + '\t' + str(varCoordNorm) + '\n')
histFile.close()
-
+
self.MPIBarrier()
def __getGlobalIndex(self, physics, iProc, iLocalVertex):
"""
Calculate the global indexing of interface nodes accross all the partitions. This does not include halo nodes.
+ This is needed because the global index of the fluid solver takes into account all the nodes, thus also those
+ in the volume mesh, not only on the interface. Here, we compute the global index of all the nodes on the
+ interface.
"""
if physics == 'fluid':
@@ -1868,254 +1879,260 @@ def __getGlobalIndex(self, physics, iProc, iLocalVertex):
globalIndex = globalStartIndex + iLocalVertex
return globalIndex
-
+
def UnsteadyFSI(self,FSI_config, FluidSolver, SolidSolver):
- """
- Run the unsteady FSI computation by synchronizing the fluid and solid solvers.
- F/s interface data are exchanged through interface mapping and interpolation (if non mathcing meshes).
- """
-
- if self.have_MPI == True:
- myid = self.comm.Get_rank()
- numberPart = self.comm.Get_size()
+ """
+ Run the unsteady FSI computation by synchronizing the fluid and solid solvers.
+ F/s interface data are exchanged through interface mapping and interpolation (if non mathcing meshes).
+ """
+
+ if self.have_MPI:
+ myid = self.comm.Get_rank()
+ numberPart = self.comm.Get_size()
else:
myid = 0
numberPart = 1
- # --- Set some general variables for the unsteady computation --- #
- deltaT = FSI_config['UNST_TIMESTEP'] # physical time step
- totTime = FSI_config['UNST_TIME'] # physical simulation time
- NbFSIIterMax = FSI_config['NB_FSI_ITER'] # maximum number of FSI iteration (for each time step)
- FSITolerance = FSI_config['FSI_TOLERANCE'] # f/s interface tolerance
- TimeIterTreshold = 0 # time iteration from which we allow the solid to deform
-
- if FSI_config['RESTART_SOL'] == 'YES':
- startTime = FSI_config['START_TIME']
- NbTimeIter = ((totTime)/deltaT)-1
- time = startTime
- TimeIter = FSI_config['RESTART_ITER']
- else:
- NbTimeIter = (totTime/deltaT)-1 # number of time iterations
- time = 0.0 # initial time
- TimeIter = 0 # initial time iteration
-
- NbTimeIter = int(NbTimeIter) # be sure that NbTimeIter is an integer
-
- varCoordNorm = 0.0 # FSI residual
- FSIConv = False # FSI convergence flag
-
- self.MPIPrint('\n**********************************')
- self.MPIPrint('* Begin unsteady FSI computation *')
- self.MPIPrint('**********************************\n')
-
- # --- Initialize the coupled solution --- #
- #If restart (DOES NOT WORK YET)
- if FSI_config['RESTART_SOL'] == 'YES':
- TimeIterTreshold = -1
- FluidSolver.setTemporalIteration(TimeIter)
- if myid == self.rootProcess:
- SolidSolver.outputDisplacements(FluidSolver.getInterRigidDispArray(), True)
- if self.have_MPI == True:
- self.comm.barrier()
- FluidSolver.setInitialMesh(True)
- if myid == self.rootProcess:
- SolidSolver.displacementPredictor(FluidSolver.getInterRigidDispArray())
- if self.have_MPI == True:
- self.comm.barrier()
- if myid == self.rootProcess:
- SolidSolver.updateSolution()
- #If no restart
- else:
- self.MPIPrint('Setting FSI initial conditions')
+ # --- Set some general variables for the unsteady computation --- #
+ deltaT = FSI_config['UNST_TIMESTEP'] # physical time step
+ totTime = FSI_config['UNST_TIME'] # physical simulation time
+ NbFSIIterMax = FSI_config['NB_FSI_ITER'] # maximum number of FSI iteration (for each time step)
+ FSITolerance = FSI_config['FSI_TOLERANCE'] # f/s interface tolerance
+ TimeIterTreshold = FSI_config['TIME_TRESHOLD'] # time iteration from which we allow the solid to deform
+ self.MPIPrint('The FSI coupling will start after {} iterations'.format(TimeIterTreshold))
+
+ if FSI_config['RESTART_SOL'] == 'YES':
+ NbTimeIter = ((totTime)/deltaT)-1
+ time = (FSI_config['RESTART_ITER']+1)*deltaT
+ TimeIter = FSI_config['RESTART_ITER']+1
+ else:
+ NbTimeIter = (totTime/deltaT)-1 # number of time iterations
+ time = 0.0 # initial time
+ TimeIter = 0 # initial time iteration
+
+ NbTimeIter = int(NbTimeIter) # be sure that NbTimeIter is an integer
+
+ varCoordNorm = 0.0 # FSI residual
+ FSIConv = False # FSI convergence flag
+
+ self.MPIPrint('\n**********************************')
+ self.MPIPrint('* Begin unsteady FSI computation *')
+ self.MPIPrint('**********************************\n')
+
+ # --- Initialize the coupled solution --- #
+ #If restart
+ if FSI_config['RESTART_SOL'] == 'YES':
+ TimeIterTreshold = -1
+ self.MPIPrint("Reading the modal amplitudes at time n-1")
if myid in self.solidSolverProcessors:
- SolidSolver.setInitialDisplacements()
+ SolidSolver.setRestart('nM1')
+ SolidSolver.setRestart('n')
self.getSolidInterfaceDisplacement(SolidSolver)
- self.interpolateSolidPositionOnFluidMesh(FSI_config)
- self.setFluidInterfaceVarCoord(FluidSolver)
- FluidSolver.SetInitialMesh() # if there is an initial deformation in the solid, it has to be communicated to the fluid solver
- self.MPIPrint('\nFSI initial conditions are set')
- self.MPIPrint('Beginning time integration\n')
-
- # --- External temporal loop --- #
- while TimeIter <= NbTimeIter:
-
- if TimeIter > TimeIterTreshold:
- NbFSIIter = NbFSIIterMax
- self.MPIPrint('\n*************** Enter Block Gauss Seidel (BGS) method for strong coupling FSI on time iteration {} ***************'.format(TimeIter))
- else:
- NbFSIIter = 1
-
- self.FSIIter = 0
- FSIConv = False
- FluidSolver.PreprocessExtIter(TimeIter) # set some parameters before temporal fluid iteration
-
- # --- Internal FSI loop --- #
- while self.FSIIter <= (NbFSIIter-1):
-
- self.MPIPrint("\n>>>> Time iteration {} / FSI iteration {} <<<<".format(TimeIter,self.FSIIter))
-
- # --- Mesh morphing step (displacements interpolation, displacements communication, and mesh morpher call) --- #
- self.interpolateSolidPositionOnFluidMesh(FSI_config)
+ self.displacementPredictor(FSI_config, SolidSolver, deltaT)
+ # We need now to update the solution because both restarter functions (solid and fluid)
+ # load the files in the solution containers, pushing back the previous solutions. We need
+ # then to push it back once more to compute the solution at the next time level
+ # Also this is required because in the fluid iteration preprocessor, if we do not update
+ # and step to the next time level, there is a flag "fsi" that will initialise the flow
+ FluidSolver.Update()
+ if myid in self.solidSolverProcessors:
+ SolidSolver.updateSolution()
+ #If no restart
+ else:
+ self.MPIPrint('Setting FSI initial conditions')
+ if myid in self.solidSolverProcessors:
+ SolidSolver.setInitialDisplacements()
+ self.getSolidInterfaceDisplacement(SolidSolver)
+ self.interpolateSolidPositionOnFluidMesh(FSI_config)
+ self.setFluidInterfaceVarCoord(FluidSolver)
+ self.MPIPrint('\nPerforming static mesh deformation (ALE) of initial mesh...\n')
+ FluidSolver.SetInitialMesh() # if there is an initial deformation in the solid, it has to be communicated to the fluid solver
+ self.MPIPrint('\nFSI initial conditions are set')
+ self.MPIPrint('Beginning time integration\n')
+
+ # --- External temporal loop --- #
+ while TimeIter <= NbTimeIter:
+
+ if TimeIter > TimeIterTreshold:
+ NbFSIIter = NbFSIIterMax
+ self.MPIPrint('\n*************** Enter Block Gauss Seidel (BGS) method for strong coupling FSI on time iteration {} ***************'.format(TimeIter))
+ else:
+ NbFSIIter = 1
+
+ self.FSIIter = 0
+ FSIConv = False
+
+ # --- Internal FSI loop --- #
+ while self.FSIIter <= (NbFSIIter-1):
+
+ self.MPIPrint("\n>>>> Time iteration {} / FSI iteration {} <<<<".format(TimeIter,self.FSIIter))
+
+ # --- Mesh morphing step (displacements interpolation, displacements communication, and mesh morpher call) --- #
+ self.interpolateSolidPositionOnFluidMesh(FSI_config)
self.MPIPrint('\nPerforming dynamic mesh deformation (ALE)...\n')
self.setFluidInterfaceVarCoord(FluidSolver)
- FluidSolver.DynamicMeshUpdate(TimeIter)
-
- # --- Fluid solver call for FSI subiteration --- #
- self.MPIPrint('\nLaunching fluid solver for one single dual-time iteration...')
- self.MPIBarrier()
- FluidSolver.ResetConvergence()
- FluidSolver.Run()
+ if self.FSIIter == 0:
+ FluidSolver.Preprocess(TimeIter) # set some parameters before temporal fluid iteration and dynamic mesh update
+ else:
+ FluidSolver.DynamicMeshUpdate(TimeIter)
+ # --- Fluid solver call for FSI subiteration --- #
+ self.MPIPrint('\nLaunching fluid solver for one single dual-time iteration...')
self.MPIBarrier()
-
- # --- Surface fluid loads interpolation and communication --- #
- self.MPIPrint('\nProcessing interface fluid loads...\n')
+ FluidSolver.ResetConvergence()
+ FluidSolver.Run()
self.MPIBarrier()
- self.getFluidInterfaceNodalForce(FSI_config, FluidSolver)
+ FluidSolver.Postprocess()
self.MPIBarrier()
- if TimeIter > TimeIterTreshold:
- self.interpolateFluidLoadsOnSolidMesh(FSI_config)
- self.setSolidInterfaceLoads(SolidSolver, FSI_config, time)
- # --- Solid solver call for FSI subiteration --- #
- self.MPIPrint('\nLaunching solid solver for a single time iteration...\n')
+ # --- Surface fluid loads interpolation and communication --- #
+ if not self.ImposedMotion:
+ self.MPIPrint('\nProcessing interface fluid loads...\n')
+ self.MPIBarrier()
+ self.getFluidInterfaceNodalForce(FSI_config, FluidSolver)
+ self.MPIBarrier()
+ if TimeIter > TimeIterTreshold:
+ if not self.ImposedMotion:
+ self.interpolateFluidLoadsOnSolidMesh(FSI_config)
+ self.setSolidInterfaceLoads(SolidSolver, FSI_config)
+
+ # --- Solid solver call for FSI subiteration --- #
+ self.MPIPrint('\nLaunching solid solver for a single time iteration...\n')
if myid in self.solidSolverProcessors:
- if FSI_config['CSD_SOLVER'] == 'NATIVE':
- SolidSolver.timeIteration(time)
- elif FSI_config['CSD_SOLVER'] == 'METAFOR' or FSI_config['CSD_SOLVER'] == 'GETDP' or FSI_config['CSD_SOLVER'] == 'TESTER':
- SolidSolver.run(time-deltaT, time)
-
- # --- Compute and monitor the FSI residual --- #
- varCoordNorm = self.computeSolidInterfaceResidual(SolidSolver)
- self.MPIPrint('\nFSI displacement norm : {}\n'.format(varCoordNorm))
- if varCoordNorm < FSITolerance:
- FSIConv = True
- break
-
- # --- Relaxe the solid position --- #
+ SolidSolver.run(time)
+
+ # --- Compute and monitor the FSI residual --- #
+ varCoordNorm = self.computeSolidInterfaceResidual(SolidSolver)
+ self.MPIPrint('\nFSI displacement norm : {}\n'.format(varCoordNorm))
+ if varCoordNorm < FSITolerance:
+ FSIConv = True
+ break
+
+ # --- Relaxe the solid position --- #
self.MPIPrint('\nProcessing interface displacements...\n')
- self.relaxSolidPosition(FSI_config)
-
- self.FSIIter += 1
- # --- End OF FSI loop --- #
+ self.relaxSolidPosition(FSI_config)
+
+ self.FSIIter += 1
+ # --- End OF FSI loop --- #
self.MPIBarrier()
- # --- Update the FSI history file --- #
- if TimeIter > TimeIterTreshold:
- self.MPIPrint('\nBGS is converged (strong coupling)')
- self.writeFSIHistory(TimeIter, time, varCoordNorm, FSIConv)
-
- # --- Update, monitor and output the fluid solution before the next time step ---#
- FluidSolver.Update()
- FluidSolver.Monitor(TimeIter)
- FluidSolver.Output(TimeIter)
-
- if TimeIter >= TimeIterTreshold:
- if myid in self.solidSolverProcessors:
- # --- Output the solid solution before thr next time step --- #
- SolidSolver.writeSolution(time, self.FSIIter, TimeIter, NbTimeIter)
-
- # --- Displacement predictor for the next time step and update of the solid solution --- #
- self.MPIPrint('\nSolid displacement prediction for next time step')
- self.displacementPredictor(FSI_config, SolidSolver, deltaT)
+ # --- Update the FSI history file --- #
+ if TimeIter > TimeIterTreshold:
+ self.MPIPrint('\nBGS is converged (strong coupling)')
+ self.writeFSIHistory(TimeIter, time, varCoordNorm, FSIConv)
+
+ # --- Update, monitor and output the fluid solution before the next time step ---#
+ FluidSolver.Update()
+ FluidSolver.Monitor(TimeIter)
+ FluidSolver.Output(TimeIter)
+
+ if TimeIter >= TimeIterTreshold:
+ if myid in self.solidSolverProcessors:
+ # --- Output the solid solution before thr next time step --- #
+ SolidSolver.writeSolution(time, TimeIter, self.FSIIter)
+
+ if TimeIter > TimeIterTreshold:
+ # --- Displacement predictor for the next time step and update of the solid solution --- #
+ self.MPIPrint('\nSolid displacement prediction for next time step')
+ self.displacementPredictor(FSI_config, SolidSolver, deltaT)
if myid in self.solidSolverProcessors:
- SolidSolver.updateSolution()
-
- TimeIter += 1
- time += deltaT
- #--- End of the temporal loop --- #
+ SolidSolver.updateSolution()
+
+ TimeIter += 1
+ time += deltaT
+ #--- End of the temporal loop --- #
self.MPIBarrier()
- self.MPIPrint('\n*************************')
- self.MPIPrint('* End FSI computation *')
- self.MPIPrint('*************************\n')
+ self.MPIPrint('\n*************************')
+ self.MPIPrint('* End FSI computation *')
+ self.MPIPrint('*************************\n')
def SteadyFSI(self, FSI_config,FluidSolver, SolidSolver):
- """
- Runs the steady FSI computation by synchronizing the fluid and solid solver with data exchange at the f/s interface.
- """
+ """
+ Runs the steady FSI computation by synchronizing the fluid and solid solver with data exchange at the f/s interface.
+ """
- if self.have_MPI == True:
- myid = self.comm.Get_rank()
- numberPart = self.comm.Get_size()
+ if self.have_MPI:
+ myid = self.comm.Get_rank()
+ numberPart = self.comm.Get_size()
else:
myid = 0
numberPart = 1
- # --- Set some general variables for the steady computation --- #
- NbIter = FSI_config['NB_EXT_ITER'] # number of fluid iteration at each FSI step
- NbFSIIterMax = FSI_config['NB_FSI_ITER'] # maximum number of FSI iteration (for each time step)
- FSITolerance = FSI_config['FSI_TOLERANCE'] # f/s interface tolerance
- varCoordNorm = 0.0
+ # --- Set some general variables for the steady computation --- #
+ NbFSIIterMax = FSI_config['NB_FSI_ITER'] # maximum number of FSI iteration (for each time step)
+ FSITolerance = FSI_config['FSI_TOLERANCE'] # f/s interface tolerance
+ varCoordNorm = 0.0
- self.MPIPrint('\n********************************')
- self.MPIPrint('* Begin steady FSI computation *')
- self.MPIPrint('********************************\n')
- self.MPIPrint('\n*************** Enter Block Gauss Seidel (BGS) method for strong coupling FSI ***************')
+ self.MPIPrint('\n********************************')
+ self.MPIPrint('* Begin steady FSI computation *')
+ self.MPIPrint('********************************\n')
+ self.MPIPrint('\n*************** Enter Block Gauss Seidel (BGS) method for strong coupling FSI ***************')
+ self.MPIPrint('Setting initial deformed mesh')
+ if myid in self.solidSolverProcessors:
+ SolidSolver.setInitialDisplacements()
self.getSolidInterfaceDisplacement(SolidSolver)
-
- # --- External FSI loop --- #
- self.FSIIter = 0
- while self.FSIIter < NbFSIIterMax:
- self.MPIPrint("\n>>>> FSI iteration {} <<<<".format(self.FSIIter))
- self.MPIPrint('\nLaunching fluid solver for a steady computation...')
- # --- Fluid solver call for FSI subiteration ---#
- Iter = 0
- FluidSolver.ResetConvergence()
- while Iter < NbIter:
- FluidSolver.PreprocessExtIter(Iter)
- FluidSolver.Run()
- StopIntegration = FluidSolver.Monitor(Iter)
- FluidSolver.Output(Iter)
- if StopIntegration:
- break;
- Iter += 1
-
- # --- Surface fluid loads interpolation and communication ---#
- self.MPIPrint('\nProcessing interface fluid loads...\n')
- self.MPIBarrier()
- self.getFluidInterfaceNodalForce(FSI_config, FluidSolver)
- self.MPIBarrier()
- self.interpolateFluidLoadsOnSolidMesh(FSI_config)
- self.setSolidInterfaceLoads(SolidSolver, FSI_config, 0.05)
-
- # --- Solid solver call for FSI subiteration --- #
- self.MPIPrint('\nLaunching solid solver for a static computation...\n')
- if myid in self.solidSolverProcessors:
- if FSI_config['CSD_SOLVER'] == 'NATIVE':
- SolidSolver.staticComputation()
- else:
- SolidSolver.run(0.0, 0.05)
- SolidSolver.writeSolution(0.0, self.FSIIter, Iter, NbIter)
-
- # --- Compute and monitor the FSI residual --- #
- varCoordNorm = self.computeSolidInterfaceResidual(SolidSolver)
- self.MPIPrint('\nFSI displacement norm : {}\n'.format(varCoordNorm))
+ self.interpolateSolidPositionOnFluidMesh(FSI_config)
+ self.setFluidInterfaceVarCoord(FluidSolver)
+ self.MPIPrint('\nFSI initial conditions are set')
+
+ # --- External FSI loop --- #
+ self.FSIIter = 0
+ while self.FSIIter < NbFSIIterMax:
+ self.MPIPrint("\n>>>> FSI iteration {} <<<<".format(self.FSIIter))
+ self.MPIPrint('\nLaunching fluid solver for a steady computation...')
+ # --- Fluid solver call for FSI subiteration ---#
+
+ FluidSolver.ResetConvergence() #This is setting to zero the convergence in the integrator, important to reset it
+ # The mesh will be deformed in the context of the preprocessor, there is no need to set the initial
+ # mesh pushing back the solution to avoid spurious velocities, as the velocity is not computed at all
+ self.MPIPrint('\nPerforming static mesh deformation...\n')
+ FluidSolver.Preprocess(0)# This will attempt to always set the initial condition, but there is a flag on the unsteady computation that will avoid it
+ FluidSolver.Run()
+ FluidSolver.Postprocess()
+ FluidSolver.Monitor(0) #This is actually not needed, it only saves the fact that the fluid solver converged innerly or reached max iterations
+ FluidSolver.Output(0)
+
+ # --- Surface fluid loads interpolation and communication ---#
+ if not self.ImposedMotion:
+ self.MPIPrint('\nProcessing interface fluid loads...\n')
+ self.MPIBarrier()
+ self.getFluidInterfaceNodalForce(FSI_config, FluidSolver)
+ self.MPIBarrier()
+ self.interpolateFluidLoadsOnSolidMesh(FSI_config)
+ self.setSolidInterfaceLoads(SolidSolver, FSI_config)
+ # --- Solid solver call for FSI subiteration --- #
+ self.MPIPrint('\nLaunching solid solver for a static computation...\n')
+ if myid in self.solidSolverProcessors:
+ SolidSolver.run(0.0)
+ SolidSolver.writeSolution(0.0, 0, self.FSIIter)
+
+ # --- Compute and monitor the FSI residual --- #
+ varCoordNorm = self.computeSolidInterfaceResidual(SolidSolver)
+ self.MPIPrint('\nFSI displacement norm : {}\n'.format(varCoordNorm))
self.writeFSIHistory(0, 0.0, varCoordNorm, False)
- if varCoordNorm < FSITolerance:
- break
+ if varCoordNorm < FSITolerance:
+ break
# --- Relaxe the solid displacement and update the solid solution --- #
self.MPIPrint('\nProcessing interface displacements...\n')
- self.relaxSolidPosition(FSI_config)
+ self.relaxSolidPosition(FSI_config)
if myid in self.solidSolverProcessors:
SolidSolver.updateSolution()
-
- # --- Mesh morphing step (displacement interpolation, displacements communication, and mesh morpher call) --- #
- self.interpolateSolidPositionOnFluidMesh(FSI_config)
- self.MPIPrint('\nPerforming static mesh deformation...\n')
- self.setFluidInterfaceVarCoord(FluidSolver)
- FluidSolver.StaticMeshUpdate()
- self.FSIIter += 1
+
+ # --- Mesh morphing step (displacement interpolation, displacements communication, and mesh morpher call) --- #
+ self.interpolateSolidPositionOnFluidMesh(FSI_config)
+ self.setFluidInterfaceVarCoord(FluidSolver)
+ self.FSIIter += 1
self.MPIBarrier()
- self.MPIPrint('\nBGS is converged (strong coupling)')
- self.MPIPrint(' ')
- self.MPIPrint('*************************')
- self.MPIPrint('* End FSI computation *')
- self.MPIPrint('*************************')
- self.MPIPrint(' ')
+ self.MPIPrint('\nBGS is converged (strong coupling)')
+ self.MPIPrint(' ')
+ self.MPIPrint('*************************')
+ self.MPIPrint('* End FSI computation *')
+ self.MPIPrint('*************************')
+ self.MPIPrint(' ')
diff --git a/SU2_PY/FSI/io/FSI_config.py b/SU2_PY/FSI_tools/FSI_config.py
similarity index 63%
rename from SU2_PY/FSI/io/FSI_config.py
rename to SU2_PY/FSI_tools/FSI_config.py
index 552a65764b68..9bc2bf9dd999 100644
--- a/SU2_PY/FSI/io/FSI_config.py
+++ b/SU2_PY/FSI_tools/FSI_config.py
@@ -2,7 +2,7 @@
## \file FSI_config.py
# \brief Python class for handling configuration file for FSI computation.
-# \author David Thomas
+# \authors Nicola Fonzi, Vittorio Cavalieri based on the work of David Thomas
# \version 7.0.8 "Blackbird"
#
# The current SU2 release has been coordinated by the
@@ -40,7 +40,7 @@
# ----------------------------------------------------------------------
import os, sys, shutil, copy
-from ..util import switch
+from FSI_tools.switch import switch
# ----------------------------------------------------------------------
# FSI Configuration Class
@@ -58,23 +58,23 @@ def __init__(self,FileName):
self.readConfig()
def __str__(self):
- tempString = str()
- for key, value in self._ConfigContent.items():
- tempString += "{} = {}\n".format(key,value)
- return tempString
+ tempString = str()
+ for key, value in self._ConfigContent.items():
+ tempString += "{} = {}\n".format(key,value)
+ return tempString
def __getitem__(self,key):
- return self._ConfigContent[key]
+ return self._ConfigContent[key]
def __setitem__(self, key, value):
- self._ConfigContent[key] = value
+ self._ConfigContent[key] = value
def readConfig(self):
input_file = open(self.ConfigFileName)
while 1:
- line = input_file.readline()
- if not line:
- break
+ line = input_file.readline()
+ if not line:
+ break
# remove line returns
line = line.strip('\r\n')
# make sure it has useful data
@@ -87,45 +87,35 @@ def readConfig(self):
for case in switch(this_param):
#integer values
- if case("NDIM") : pass
- #if case("MESH_DEF_LIN_ITER") : pass
- #if case("MESH_DEF_NONLIN_ITER") : pass
- if case("RESTART_ITER") : pass
- if case("NB_EXT_ITER") : pass
- if case("NB_FSI_ITER") :
- self._ConfigContent[this_param] = int(this_value)
- break
-
- #float values
+ if case("NDIM") : pass
+ if case("RESTART_ITER") : pass
+ if case("TIME_TRESHOLD") : pass
+ if case("NB_FSI_ITER") :
+ self._ConfigContent[this_param] = int(this_value)
+ break
+
+ #float values
if case("RBF_RADIUS") : pass
- if case("AITKEN_PARAM") : pass
- if case("START_TIME") : pass
- if case("UNST_TIMESTEP") : pass
- if case("UNST_TIME") : pass
- if case("FSI_TOLERANCE") :
- self._ConfigContent[this_param] = float(this_value)
- break
-
- #string values
- if case("CFD_CONFIG_FILE_NAME") : pass
- if case("CSD_SOLVER") : pass
- if case("CSD_CONFIG_FILE_NAME") : pass
- if case("RESTART_SOL") : pass
- if case("MATCHING_MESH") : pass
+ if case("AITKEN_PARAM") : pass
+ if case("UNST_TIMESTEP") : pass
+ if case("UNST_TIME") : pass
+ if case("FSI_TOLERANCE") :
+ self._ConfigContent[this_param] = float(this_value)
+ break
+
+ #string values
+ if case("CFD_CONFIG_FILE_NAME") : pass
+ if case("CSD_SOLVER") : pass
+ if case("CSD_CONFIG_FILE_NAME") : pass
+ if case("RESTART_SOL") : pass
+ if case("MATCHING_MESH") : pass
if case("MESH_INTERP_METHOD") : pass
- if case("DISP_PRED") : pass
- if case("AITKEN_RELAX") : pass
- if case("TIME_MARCHING") : pass
- if case("INTERNAL_FLOW") :
- #if case("MESH_DEF_METHOD") : pass
- self._ConfigContent[this_param] = this_value
- break
-
- if case():
- print(this_param + " is an invalid option !")
- break
- #end for
-
-
-
- #def dump()
+ if case("DISP_PRED") : pass
+ if case("AITKEN_RELAX") : pass
+ if case("TIME_MARCHING") :
+ self._ConfigContent[this_param] = this_value
+ break
+
+ if case():
+ print(this_param + " is an invalid option !")
+ break
diff --git a/SU2_PY/FSI_tools/__init__.py b/SU2_PY/FSI_tools/__init__.py
new file mode 100644
index 000000000000..e843f12bbe49
--- /dev/null
+++ b/SU2_PY/FSI_tools/__init__.py
@@ -0,0 +1,3 @@
+from FSI_tools.FSIInterface import Interface
+from FSI_tools.switch import switch
+from FSI_tools.FSI_config import FSIConfig
diff --git a/SU2_PY/FSI/util/switch.py b/SU2_PY/FSI_tools/switch.py
similarity index 96%
rename from SU2_PY/FSI/util/switch.py
rename to SU2_PY/FSI_tools/switch.py
index 08fef0f24463..b42eaf6e9ddd 100644
--- a/SU2_PY/FSI/util/switch.py
+++ b/SU2_PY/FSI_tools/switch.py
@@ -1,52 +1,52 @@
-# -------------------------------------------------------------------
-# Switch Class
-# -------------------------------------------------------------------
-# source: Brian Beck, PSF License, ActiveState Code
-# http://code.activestate.com/recipes/410692/
-
-class switch(object):
- """ Readable switch construction
-
- Example:
-
- c = 'z'
- for case in switch(c):
- if case('a'): pass # only necessary if the rest of the suite is empty
- if case('b'): pass
- # ...
- if case('y'): pass
- if case('z'):
- print("c is lowercase!")
- break
- if case('A'): pass
- # ...
- if case('Z'):
- print("c is uppercase!")
- break
- if case(): # default
- print("I dunno what c was!")
-
- source: Brian Beck, PSF License, ActiveState Code
- http://code.activestate.com/recipes/410692/
- """
-
- def __init__(self, value):
- self.value = value
- self.fall = False
-
- def __iter__(self):
- """Return the match method once, then stop"""
- yield self.match
- raise StopIteration
-
- def match(self, *args):
- """Indicate whether or not to enter a case suite"""
- if self.fall or not args:
- return True
- elif self.value in args:
- self.fall = True
- return True
- else:
- return False
-
-#: class switch()
+# -------------------------------------------------------------------
+# Switch Class
+# -------------------------------------------------------------------
+# source: Brian Beck, PSF License, ActiveState Code
+# http://code.activestate.com/recipes/410692/
+
+class switch(object):
+ """ Readable switch construction
+
+ Example:
+
+ c = 'z'
+ for case in switch(c):
+ if case('a'): pass # only necessary if the rest of the suite is empty
+ if case('b'): pass
+ # ...
+ if case('y'): pass
+ if case('z'):
+ print("c is lowercase!")
+ break
+ if case('A'): pass
+ # ...
+ if case('Z'):
+ print("c is uppercase!")
+ break
+ if case(): # default
+ print("I dunno what c was!")
+
+ source: Brian Beck, PSF License, ActiveState Code
+ http://code.activestate.com/recipes/410692/
+ """
+
+ def __init__(self, value):
+ self.value = value
+ self.fall = False
+
+ def __iter__(self):
+ """Return the match method once, then stop"""
+ yield self.match
+ raise StopIteration
+
+ def match(self, *args):
+ """Indicate whether or not to enter a case suite"""
+ if self.fall or not args:
+ return True
+ elif self.value in args:
+ self.fall = True
+ return True
+ else:
+ return False
+
+#: class switch()
diff --git a/SU2_PY/SU2_Nastran/__init__.py b/SU2_PY/SU2_Nastran/__init__.py
new file mode 100644
index 000000000000..73eb7e0b879b
--- /dev/null
+++ b/SU2_PY/SU2_Nastran/__init__.py
@@ -0,0 +1 @@
+from SU2_Nastran import pysu2_nastran
diff --git a/SU2_PY/SU2_Nastran/pysu2_nastran.py b/SU2_PY/SU2_Nastran/pysu2_nastran.py
new file mode 100644
index 000000000000..2d6185a2149b
--- /dev/null
+++ b/SU2_PY/SU2_Nastran/pysu2_nastran.py
@@ -0,0 +1,909 @@
+#!/usr/bin/env python
+
+## \file pysu2_nastran.py
+# \brief Structural solver using Nastran models
+# \authors Nicola Fonzi, Vittorio Cavalieri, based on the work of David Thomas
+# \version 7.0.8 "Blackbird"
+#
+# SU2 Project Website: https://su2code.github.io
+#
+# The SU2 Project is maintained by the SU2 Foundation
+# (http://su2foundation.org)
+#
+# Copyright 2012-2020, 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 .
+
+# ----------------------------------------------------------------------
+# Imports
+# ----------------------------------------------------------------------
+
+import os, sys, shutil, copy
+import numpy as np
+import scipy as sp
+import scipy.linalg as linalg
+from math import *
+from FSI_tools.switch import switch
+
+# ----------------------------------------------------------------------
+# Config class
+# ----------------------------------------------------------------------
+
+class RefSystem:
+
+ def __init__(self):
+ self.CID = 0
+ self.RID = 0
+ self.Origin = np.array([[0.],[0.],[0.]])
+ self.Rot = np.array([[0.,0.,0.],[0.,0.,0.],[0.,0.,0.]])
+
+ def SetOrigin(self,A):
+ AX , AY , AZ = A
+ self.Origin[0] = AX
+ self.Origin[1] = AY
+ self.Origin[2] = AZ
+
+ def SetRotMatrix(self,x,y,z):
+ self.Rot = np.array([[x[0],y[0],z[0]],[x[1],y[1],z[1]],[x[2],y[2],z[2]]])
+
+ def SetCID(self,CID):
+ self.CID = CID
+
+ def SetRID(self,RID):
+ self.RID = RID
+
+ def GetOrigin(self):
+ return self.Origin
+
+ def GetRotMatrix(self):
+ return self.Rot
+
+ def GetRID(self):
+ return self.RID
+
+ def GetCID(self):
+ return self.CID
+
+class Point:
+ """
+ Class containing data regarding all the structural nodes.
+ Coord0: Coordinates at the initial time iteration.
+ Coord: Coordinates at the current time iteration.
+ Coord_n: Coordinates at the previous time iteration.
+ Vel: Velocity at the current time iteration.
+ Vel_n: Velocity at the previous time iteration.
+ Force: Nodal force provided by the aerodynamics.
+ ID: ID of the node.
+ CP: Coordinate system definition of the position.
+ CD: Coordinate system definition of the output coming from Nastran.
+ """
+
+ def __init__(self):
+ self.Coord0 = np.zeros((3,1))
+ self.Coord = np.zeros((3,1))
+ self.Coord_n = np.zeros((3,1))
+ self.Vel = np.zeros((3,1))
+ self.Vel_n = np.zeros((3,1))
+ self.Force = np.zeros((3,1))
+ self.ID = 0
+ self.CP = 0
+ self.CD = 0
+
+ def GetCoord0(self):
+ return self.Coord0
+
+ def GetCoord(self):
+ return self.Coord
+
+ def GetCoord_n(self):
+ return self.Coord_n
+
+ def GetVel(self):
+ return self.Vel
+
+ def GetVel_n(self):
+ return self.Vel_n
+
+ def GetForce(self):
+ return self.Force
+
+ def GetID(self):
+ return self.ID
+
+ def GetCP(self):
+ return self.CP
+
+ def GetCD(self):
+ return self.CD
+
+ def SetCoord0(self, val_Coord):
+ x, y, z = val_Coord
+ self.Coord0[0] = x
+ self.Coord0[1] = y
+ self.Coord0[2] = z
+
+ def SetCoord(self, val_Coord):
+ x, y, z = val_Coord
+ self.Coord[0] = x
+ self.Coord[1] = y
+ self.Coord[2] = z
+
+ def SetCoord_n(self, val_Coord):
+ x, y, z = val_Coord
+ self.Coord_n[0] = x
+ self.Coord_n[1] = y
+ self.Coord_n[2] = z
+
+ def SetVel(self, val_Vel):
+ vx, vy, vz = val_Vel
+ self.Vel[0] = vx
+ self.Vel[1] = vy
+ self.Vel[2] = vz
+
+ def SetVel_n(self, val_Vel):
+ vx, vy, vz = val_Vel
+ self.Vel_n[0] = vx
+ self.Vel_n[1] = vy
+ self.Vel_n[2] = vz
+
+ def SetForce(self, val_Force):
+ fx, fy, fz = val_Force
+ self.Force[0] = fx
+ self.Force[1] = fy
+ self.Force[2] = fz
+
+ def SetID(self, ID):
+ self.ID = ID
+
+ def SetCP(self,CP):
+ self.CP = CP
+
+ def SetCD(self,CD):
+ self.CD = CD
+
+ def updateCoordVel(self):
+ self.Coord_n = np.copy(self.Coord)
+ self.Vel_n = np.copy(self.Vel)
+
+class Solver:
+ """
+ Structural solver main class.
+ It contains all the required methods for the coupling with SU2.
+ """
+
+ def __init__(self, config_fileName, ImposedMotion):
+ """
+ Constructor of the structural solver class.
+ """
+
+ self.Config_file = config_fileName
+ self.Config = {}
+
+ print("\n------------------------------ Configuring the structural tester solver for FSI simulation ------------------------------")
+ self.__readConfig()
+
+ self.Mesh_file = self.Config['MESH_FILE']
+ self.Punch_file = self.Config['PUNCH_FILE']
+ self.FSI_marker = self.Config['MOVING_MARKER']
+ self.Unsteady = (self.Config['TIME_MARCHING']=="YES")
+ self.ImposedMotion = ImposedMotion
+ if self.Unsteady:
+ print('Dynamic computation.')
+ self.nDof = self.Config['NMODES']
+ print("Reading number of modes from file")
+
+
+ # Structural properties
+ print("Reading the modal and stiffnes matrix from file")
+ self.ModalDamping = self.Config['MODAL_DAMPING']
+ if self.ModalDamping == 0:
+ print("The structural model is undamped")
+ else:
+ print("Assuming {}% of modal damping".format(self.ModalDamping*100))
+
+ self.deltaT = self.Config['DELTA_T']
+ self.rhoAlphaGen = self.Config['RHO']
+
+ self.nElem = int()
+ self.nPoint = int()
+ self.nMarker = int()
+ self.nRefSys = int()
+ self.node = []
+ self.markers = {}
+ self.refsystems = []
+
+ print("\n------------------------------ Reading the mesh ------------------------------")
+ self.__readNastranMesh()
+
+ print("\n------------------------------ Creating the structural model ------------------------------")
+ self.__setStructuralMatrices()
+
+ print("\n------------------------------ Setting the integration parameters ------------------------------")
+ self.__setIntegrationParameters()
+ self.__setInitialConditions()
+
+ # Prepare the output file
+ if self.Config["RESTART_SOL"]=="NO":
+ histFile = open('StructHistoryModal.dat', "w")
+ header = 'Time\t' + 'Time Iteration\t' + 'FSI Iteration\t'
+ for imode in range(self.nDof):
+ header = header + 'q' + str(imode+1) + '\t' + 'qdot' + str(imode+1) + '\t' + 'qddot' + str(imode+1) + '\t'
+ header = header + '\n'
+ histFile.write(header)
+ histFile.close()
+
+ def __readConfig(self):
+ """
+ This methods obtains the configuration options from the structural solver input
+ file.
+ """
+
+ with open(self.Config_file) as configfile:
+ while 1:
+ line = configfile.readline()
+ if not line:
+ break
+
+ # remove line returns
+ line = line.strip('\r\n')
+ # make sure it has useful data
+ if (not "=" in line) or (line[0] == '%'):
+ continue
+ # split across equal sign
+ line = line.split("=",1)
+ this_param = line[0].strip()
+ this_value = line[1].strip()
+
+ for case in switch(this_param):
+ #integer values
+ if case("NMODES") : pass
+ if case("IMPOSED_MODE") : pass
+ if case("RESTART_ITER") :
+ self.Config[this_param] = int(this_value)
+ break
+
+ #float values
+ if case("DELTA_T") : pass
+ if case("MODAL_DAMPING") : pass
+ if case("RHO") :
+ self.Config[this_param] = float(this_value)
+ break
+
+ #string values
+ if case("TIME_MARCHING") : pass
+ if case("MESH_FILE") : pass
+ if case("PUNCH_FILE") : pass
+ if case("RESTART_SOL") : pass
+ if case("IMPOSED_DISP") : pass
+ if case("IMPOSED_VEL") : pass
+ if case("IMPOSED_ACC") : pass
+ if case("MOVING_MARKER") :
+ self.Config[this_param] = this_value
+ break
+
+ #lists values
+ if case("INITIAL_MODES"):
+ self.Config[this_param] = eval(this_value)
+ break
+
+ if case():
+ print(this_param + " is an invalid option !")
+ break
+
+
+
+ def __readNastranMesh(self):
+ """
+ This method reads the nastran 3D mesh.
+ """
+
+ def nastran_float(s):
+ if s.find('E') == -1:
+ s = s.replace('-','e-')
+ s = s.replace('+','e+')
+ if s[0] == 'e':
+ s = s[1:]
+ return float(s)
+
+ self.nMarker = 1
+ self.nPoint = 0
+ self.nRefSys = 0
+
+ with open(self.Mesh_file,'r') as meshfile:
+ print('Opened mesh file ' + self.Mesh_file + '.')
+ while 1:
+ line = meshfile.readline()
+ if not line:
+ break
+
+ pos = line.find('GRID')
+ if pos == 30:
+ line = line.strip('\r\n')
+ self.node.append(Point())
+ line = line[30:]
+ ID = int(line[8:16])
+ CP = int(line[16:24])
+ x = nastran_float(line[24:32])
+ y = nastran_float(line[32:40])
+ z = nastran_float(line[40:48])
+ if CP != 0:
+ for iRefSys in range(self.nRefSys):
+ if self.refsystems[iRefSys].GetCID()==CP:
+ break
+ if self.refsystems[iRefSys].GetCID()!=CP:
+ sys.exit('Definition reference {} system not found'.format(CP))
+ DeltaPos = self.refsystems[iRefSys].GetOrigin()
+ RotatedPos = self.refsystems[iRefSys].GetRotMatrix().dot(np.array([[x],[y],[z]]))
+ x = RotatedPos[0]+DeltaPos[0]
+ y = RotatedPos[1]+DeltaPos[1]
+ z = RotatedPos[2]+DeltaPos[2]
+ CD = int(line[48:56])
+ self.node[self.nPoint].SetCoord((x,y,z))
+ self.node[self.nPoint].SetID(ID)
+ self.node[self.nPoint].SetCP(CP)
+ self.node[self.nPoint].SetCD(CD)
+ self.node[self.nPoint].SetCoord0((x,y,z))
+ self.node[self.nPoint].SetCoord_n((x,y,z))
+ self.nPoint = self.nPoint+1
+ continue
+
+ pos = line.find('CORD2R')
+ if pos == 30:
+ line = line.strip('\r\n')
+ self.refsystems.append(RefSystem())
+ line = line[30:]
+ CID = int(line[8:16])
+ self.refsystems[self.nRefSys].SetCID(CID)
+ RID = int(line[16:24])
+ if RID!=0:
+ print('ERROR: Reference system {} must be defined with respect to global reference system'.format(CID))
+ sys.exit()
+ self.refsystems[self.nRefSys].SetRID(RID)
+ AX = nastran_float(line[24:32])
+ AY = nastran_float(line[32:40])
+ AZ = nastran_float(line[40:48])
+ BX = nastran_float(line[48:56])
+ BY = nastran_float(line[56:64])
+ BZ = nastran_float(line[64:72])
+ z_direction = np.array([BX-AX,BY-AY,BZ-AZ])
+ z_direction = z_direction/linalg.norm(z_direction)
+ line = meshfile.readline()
+ line = line.strip('\r\n')
+ line = line[30:]
+ CX = nastran_float(line[8:16])
+ CY = nastran_float(line[16:24])
+ CZ = nastran_float(line[24:32])
+ y_direction = np.cross(z_direction,[CX-AX,CY-AY,CZ-AZ])
+ y_direction = y_direction/linalg.norm(y_direction)
+ x_direction = np.cross(y_direction,z_direction)
+ x_direction = x_direction/linalg.norm(x_direction)
+ self.refsystems[self.nRefSys].SetRotMatrix(x_direction,y_direction,z_direction)
+ self.refsystems[self.nRefSys].SetOrigin((AX,AY,AZ))
+ self.nRefSys = self.nRefSys+1
+ continue
+
+ pos = line.find("SET1")
+ markerTag = self.FSI_marker
+ if pos == 30:
+ self.markers[markerTag] = []
+ line = line.strip('\r\n')
+ line = line[46:]
+ line = line.split()
+ existValue = True
+ while existValue:
+ if line[0] == "+":
+ line = meshfile.readline()
+ line = line.strip('\r\n')
+ line = line[37:]
+ line = line.split()
+ ID = int(line.pop(0))
+ for iPoint in range(self.nPoint):
+ if self.node[iPoint].GetID() == ID:
+ break
+ self.markers[self.FSI_marker].append(iPoint)
+ existValue = len(line)>=1
+ continue
+
+ self.markers[self.FSI_marker].sort()
+ print("Number of elements: {}".format(self.nElem))
+ print("Number of point: {}".format(self.nPoint))
+ print("Number of markers: {}".format(self.nMarker))
+ print("Number of reference systems: {}".format(self.nRefSys))
+ if len(self.markers) > 0:
+ print("Moving marker(s):")
+ for mark in self.markers.keys():
+ print(mark)
+
+ def __setStructuralMatrices(self):
+ """
+ This method reads the punch file and obtains the modal shapes and modal stiffnesses.
+ """
+
+ self.M = np.zeros((self.nDof, self.nDof))
+ self.K = np.zeros((self.nDof, self.nDof))
+ self.C = np.zeros((self.nDof, self.nDof))
+
+ self.q = np.zeros((self.nDof, 1))
+ self.qdot = np.zeros((self.nDof, 1))
+ self.qddot = np.zeros((self.nDof, 1))
+ self.a = np.zeros((self.nDof, 1))
+
+ self.q_n = np.zeros((self.nDof, 1))
+ self.qdot_n = np.zeros((self.nDof, 1))
+ self.qddot_n = np.zeros((self.nDof, 1))
+ self.a_n = np.zeros((self.nDof, 1))
+
+ self.F = np.zeros((self.nDof, 1))
+
+ self.Ux = np.zeros((self.nPoint,self.nDof))
+ self.Uy = np.zeros((self.nPoint,self.nDof))
+ self.Uz = np.zeros((self.nPoint,self.nDof))
+
+ with open(self.Punch_file,'r') as punchfile:
+ print('Opened punch file ' + self.Punch_file + '.')
+ while 1:
+ line = punchfile.readline()
+ if not line:
+ break
+
+ pos = line.find('MODE ')
+ if pos != -1:
+ line = line.strip('\r\n').split()
+ n = int(line[5])
+ imode = n-1
+ k_i = float(line[2])
+ self.M[imode][imode] = 1
+ self.K[imode][imode] = k_i
+ w_i = sqrt(k_i)
+ self.C[imode][imode] = 2 * self.ModalDamping * w_i
+ iPoint = 0
+ for indexIter in range(self.nPoint):
+ line = punchfile.readline()
+ line = line.strip('\r\n').split()
+ if line[1]=='G':
+ ux = float(line[2])
+ uy = float(line[3])
+ uz = float(line[4])
+ if self.node[iPoint].GetCD()!=0:
+ for iRefSys in range(self.nRefSys):
+ if self.refsystems[iRefSys].GetCID()==self.node[iPoint].GetCD():
+ break
+ if self.refsystems[iRefSys].GetCID()!=self.node[iPoint].GetCD():
+ sys.exit('Output reference {} system not found'.format(self.node[iPoint].GetCD()))
+ RotatedOutput = self.refsystems[iRefSys].GetRotMatrix().dot(np.array([[ux],[uy],[uz]]))
+ ux = RotatedOutput[0]
+ uy = RotatedOutput[1]
+ uz = RotatedOutput[2]
+ self.Ux[iPoint][imode] = ux
+ self.Uy[iPoint][imode] = uy
+ self.Uz[iPoint][imode] = uz
+ iPoint = iPoint + 1
+ line = punchfile.readline()
+ if line[1]=='S':
+ line = punchfile.readline()
+
+ if n == self.nDof:
+ break
+
+ self.__setNonDiagonalStructuralMatrices()
+
+ self.UxT = self.Ux.transpose()
+ self.UyT = self.Uy.transpose()
+ self.UzT = self.Uz.transpose()
+
+ if n= eps:
+ St = self.__TangentOperator()
+ Deltaq = -1*(linalg.solve(St,res))
+ self.q += Deltaq
+ self.qdot += self.gammaPrime*Deltaq
+ self.qddot += self.betaPrime*Deltaq
+ res = self.__ComputeResidual()
+
+ self.a += (1-self.alpha_f)/(1-self.alpha_m)*self.qddot
+ else:
+ self.q[self.Config["IMPOSED_MODE"]] = eval(self.Config["IMPOSED_DISP"])
+ self.qdot[self.Config["IMPOSED_MODE"]] = eval(self.Config["IMPOSED_VEL"])
+ self.qddot[self.Config["IMPOSED_MODE"]] = eval(self.Config["IMPOSED_ACC"])
+ self.a = np.copy(self.qddot)
+
+
+ def __SetLoads(self):
+ """
+ This method uses the nodal forces and the mode shapes to obtain the modal forces.
+ """
+ makerID = list(self.markers.keys())
+ makerID = makerID[0]
+ nodeList = self.markers[makerID]
+ FX = np.zeros((self.nPoint, 1))
+ FY = np.zeros((self.nPoint, 1))
+ FZ = np.zeros((self.nPoint, 1))
+ for iPoint in nodeList:
+ Force = self.node[iPoint].GetForce()
+ FX[iPoint] = float(Force[0])
+ FY[iPoint] = float(Force[1])
+ FZ[iPoint] = float(Force[2])
+ self.F = self.UxT.dot(FX) + self.UyT.dot(FY) + self.UzT.dot(FZ)
+
+ def __ComputeResidual(self):
+ """
+ This method computes the residual for integration.
+ """
+
+ res = self.M.dot(self.qddot) + self.C.dot(self.qdot) + self.K.dot(self.q) - self.F
+
+ return res
+
+ def __TangentOperator(self):
+ """
+ This method computes the tangent operator for solution.
+ """
+
+ # The problem is linear, so the tangent operator is straightforward.
+ St = self.betaPrime*self.M + self.gammaPrime*self.C + self.K
+
+ return St
+
+ def exit(self):
+ """
+ This method cleanly exits the structural solver.
+ """
+
+ print("\n**************** Exiting the structural tester solver ****************")
+
+ def run(self,time):
+ """
+ This method is the main function for advancing the solution of one time step.
+ """
+ self.__temporalIteration(time)
+ header = 'Time\t'
+ for imode in range(min([self.nDof,5])):
+ header = header + 'q' + str(imode+1) + '\t' + 'qdot' + str(imode+1) + '\t' + 'qddot' + str(imode+1) + '\t'
+ header = header + '\n'
+ print(header)
+ line = '{:6.4f}'.format(time) + '\t'
+ for imode in range(min([self.nDof,5])):
+ line = line + '{:6.4f}'.format(float(self.q[imode])) + '\t' + '{:6.4f}'.format(float(self.qdot[imode])) + '\t' + '{:6.4f}'.format(float(self.qddot[imode])) + '\t'
+ line = line + '\n'
+ print(line)
+ self.__computeInterfacePosVel(False)
+
+ def setInitialDisplacements(self):
+ """
+ This method provides public access to the method __computeInterfacePosVel and
+ sets velocities for previous time steps.
+ """
+
+ self.__computeInterfacePosVel(True)
+
+ def setRestart(self, timeIter):
+ if timeIter == 'nM1':
+ #read the Structhistory to obtain the mode amplitudes
+ with open('StructHistoryModal.dat','r') as file:
+ print('Opened history file ' + 'StructHistoryModal.dat' + '.')
+ line = file.readline()
+ while 1:
+ line = file.readline()
+ if not line:
+ print("The restart iteration was not found in the structural history")
+ break
+ line = line.strip('\r\n').split()
+ if int(line[1])==(self.Config["RESTART_ITER"]-1):
+ break
+ index = 0
+ for index_mode in range(self.nDof):
+ self.q[index_mode] = float(line[index+3])
+ self.qdot[index_mode] = float(line[index+4])
+ self.qddot[index_mode] = float(line[index+5])
+ index += 3
+ del index
+ #push back the mode amplitudes velocities and accelerations
+ self.__computeInterfacePosVel(True)
+ self.q_n = np.copy(self.q)
+ self.qdot_n = np.copy(self.qdot)
+ self.qddot_n = np.copy(self.qddot)
+ self.a_n = np.copy(self.a)
+ if timeIter == 'n':
+ #read the Structhistory to obtain the modes
+ with open('StructHistoryModal.dat','r') as file:
+ print('Opened history file ' + 'StructHistoryModal.dat' + '.')
+ line = file.readline()
+ while 1:
+ line = file.readline()
+ if not line:
+ print("The restart iteration was not found in the structural history")
+ break
+ line = line.strip('\r\n').split()
+ if int(line[1])==self.Config["RESTART_ITER"]:
+ break
+ index = 0
+ for index_mode in range(self.nDof):
+ self.q[index_mode] = float(line[index+3])
+ self.qdot[index_mode] = float(line[index+4])
+ self.qddot[index_mode] = float(line[index+5])
+ index += 3
+ del index
+ self.__computeInterfacePosVel(False)
+
+ def writeSolution(self, time, timeIter, FSIIter):
+ """
+ This method is the main function for output. It writes the file StructHistoryModal.dat
+ """
+
+ # Modal History
+ histFile = open('StructHistoryModal.dat', "a")
+ line = str(time) + '\t' + str(timeIter) + '\t' + str(FSIIter) + '\t'
+ for imode in range(self.nDof):
+ line = line + str(float(self.q[imode])) + '\t' + str(float(self.qdot[imode])) + '\t' + str(float(self.qddot[imode])) + '\t'
+ line = line + '\n'
+ histFile.write(line)
+ histFile.close()
+
+ def updateSolution(self):
+ """
+ This method updates the solution.
+ """
+
+ self.q_n = np.copy(self.q)
+ self.qdot_n = np.copy(self.qdot)
+ self.qddot_n = np.copy(self.qddot)
+ self.a_n = np.copy(self.a)
+ self.__reset(self.q)
+ self.__reset(self.qdot)
+ self.__reset(self.qddot)
+ self.__reset(self.a)
+
+ makerID = list(self.markers.keys())
+ makerID = makerID[0]
+ nodeList = self.markers[makerID]
+
+ for iPoint in nodeList:
+ self.node[iPoint].updateCoordVel()
+
+
+ def applyload(self, iVertex, fx, fy, fz):
+ """
+ This method can be accessed from outside to set the nodal forces.
+ """
+
+ makerID = list(self.markers.keys())
+ makerID = makerID[0]
+ iPoint = self.getInterfaceNodeGlobalIndex(makerID, iVertex)
+ self.node[iPoint].SetForce((fx,fy,fz))
+
+ def getFSIMarkerID(self):
+ """
+ This method provides the ID of the interface marker
+ """
+ L = list(self.markers)
+ return L[0]
+
+ def getNumberOfSolidInterfaceNodes(self, markerID):
+
+ return len(self.markers[markerID])
+
+ def getInterfaceNodeGlobalIndex(self, markerID, iVertex):
+
+ return self.markers[markerID][iVertex]
+
+ def getInterfaceNodePos(self, markerID, iVertex):
+
+ iPoint = self.markers[markerID][iVertex]
+ Coord = self.node[iPoint].GetCoord()
+ return Coord
+
+ def getInterfaceNodeDisp(self, markerID, iVertex):
+
+ iPoint = self.markers[markerID][iVertex]
+ Coord = self.node[iPoint].GetCoord()
+ Coord0 = self.node[iPoint].GetCoord0()
+ return (Coord-Coord0)
+
+ def getInterfaceNodeVel(self, markerID, iVertex):
+
+ iPoint = self.markers[markerID][iVertex]
+ Vel = self.node[iPoint].GetVel()
+ return Vel
+
+ def getInterfaceNodeVelNm1(self, markerID, iVertex):
+
+ iPoint = self.markers[markerID][iVertex]
+ Vel = self.node[iPoint].GetVel_n()
+ return Vel
diff --git a/SU2_PY/fsi_computation.py b/SU2_PY/fsi_computation.py
old mode 100755
new mode 100644
index a2e1350ea69d..c4b8a0e277fc
--- a/SU2_PY/fsi_computation.py
+++ b/SU2_PY/fsi_computation.py
@@ -2,12 +2,12 @@
## \file fsi_computation.py
# \brief Python wrapper code for FSI computation by coupling a third-party structural solver to SU2.
-# \author David Thomas
+# \authors Nicola Fonzi, Vittorio Cavalieri based on the work of David Thomas
# \version 7.0.8 "Blackbird"
#
# SU2 Project Website: https://su2code.github.io
-#
-# The SU2 Project is maintained by the SU2 Foundation
+#
+# The SU2 Project is maintained by the SU2 Foundation
# (http://su2foundation.org)
#
# Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md)
@@ -16,7 +16,7 @@
# 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
@@ -37,15 +37,12 @@
from math import * # use mathematical expressions
from optparse import OptionParser # use a parser for configuration
-import SU2 # imports SU2 python tools
-import FSI # imports FSI python tools
-
-
# imports the CFD (SU2) module for FSI computation
import pysu2
+import FSI_tools as FSI # imports FSI python tools
# -------------------------------------------------------------------
-# Main
+# Main
# -------------------------------------------------------------------
def main():
@@ -76,17 +73,20 @@ def main():
# --- Set the working directory --- #
if myid == rootProcess:
if os.getcwd() not in sys.path:
- sys.path.append(os.getcwd())
- print("Setting working directory : {}".format(os.getcwd()))
+ sys.path.append(os.getcwd())
+ print("Setting working directory : {}".format(os.getcwd()))
else:
- print("Working directory is set to {}".format(os.getcwd()))
+ print("Working directory is set to {}".format(os.getcwd()))
+
+ if have_MPI:
+ comm.barrier()
# starts timer
start = timer.time()
confFile = str(options.filename)
- FSI_config = FSI.io.FSIConfig(confFile) # FSI configuration file
+ FSI_config = FSI.FSIConfig(confFile) # FSI configuration file
CFD_ConFile = FSI_config['CFD_CONFIG_FILE_NAME'] # CFD configuration file
CSD_ConFile = FSI_config['CSD_CONFIG_FILE_NAME'] # CSD configuration file
@@ -99,9 +99,9 @@ def main():
if myid == rootProcess:
print('\n***************************** Initializing fluid solver *****************************')
try:
- FluidSolver = pysu2.CFluidDriver(CFD_ConFile, 1, FSI_config['NDIM'], comm)
+ FluidSolver = pysu2.CSinglezoneDriver(CFD_ConFile, 1, comm)
except TypeError as exception:
- print('A TypeError occured in pysu2.CSingleZoneDriver : ',exception)
+ print('A TypeError occured in pysu2.CSinglezoneDriver : ',exception)
if have_MPI:
print('ERROR : You are trying to initialize MPI with a serial build of the wrapper. Please, remove the --parallel option that is incompatible with a serial build.')
else:
@@ -110,21 +110,18 @@ def main():
if have_MPI:
comm.barrier()
-
+
# --- Initialize the solid solver --- # (!! for now we are using only serial solid solvers)
if myid == rootProcess:
print('\n***************************** Initializing solid solver *****************************')
- if CSD_Solver == 'METAFOR':
- from MetaforSolver import MtfSolver
- SolidSolver = MtfSolver(CSD_ConFile)
- elif CSD_Solver == 'NATIVE':
- import NativeSolid
- SolidSolver = NativeSolid.NativeSolidSolver(CSD_ConFile, True)
- elif CSD_Solver == 'GETDP':
- import GetDPSolver
- SolidSolver = GetDPSolver.GetDPSolver(CSD_ConFile, True)
- elif CSD_Solver == 'TESTER':
- SolidSolver = FSI.PitchPlungeAirfoilStructuralTester.Solver(CSD_ConFile)
+ if CSD_Solver == 'AEROELASTIC':
+ from SU2_Nastran import pysu2_nastran
+ SolidSolver = pysu2_nastran.Solver(CSD_ConFile,False)
+ elif CSD_Solver == 'IMPOSED':
+ from SU2_Nastran import pysu2_nastran
+ SolidSolver = pysu2_nastran.Solver(CSD_ConFile,True)
+ else:
+ print("\n Invalid solid solver option")
else:
SolidSolver = None
@@ -137,7 +134,7 @@ def main():
if have_MPI:
comm.barrier()
FSIInterface = FSI.Interface(FSI_config, FluidSolver, SolidSolver, have_MPI)
-
+
if myid == rootProcess:
print('\n***************************** Connect fluid and solid solvers *****************************')
if have_MPI:
@@ -149,7 +146,7 @@ def main():
if have_MPI:
comm.barrier()
FSIInterface.interfaceMapping(FluidSolver, SolidSolver, FSI_config)
-
+
if have_MPI:
comm.barrier()
@@ -168,7 +165,6 @@ def main():
print('A KeyboardInterrupt occured in FSIInterface.UnsteadyFSI : ',exception)
else:
try:
- NbExtIter = FSI_config['NB_EXT_ITER']
FSIInterface.SteadyFSI(FSI_config, FluidSolver, SolidSolver)
except NameError as exception:
if myid == rootProcess:
@@ -179,7 +175,7 @@ def main():
except KeyboardInterrupt as exception :
if myid == rootProcess:
print('A KeyboardInterrupt occured in FSIInterface.SteadyFSI : ',exception)
-
+
if have_MPI:
comm.barrier()
@@ -194,7 +190,7 @@ def main():
# stops timer
stop = timer.time()
elapsedTime = stop-start
-
+
if myid == rootProcess:
print("\n Computation successfully performed in {} seconds.".format(elapsedTime))
diff --git a/SU2_PY/meson.build b/SU2_PY/meson.build
index 319fad083049..fa11084efc76 100644
--- a/SU2_PY/meson.build
+++ b/SU2_PY/meson.build
@@ -13,7 +13,7 @@ install_data(['continuous_adjoint.py',
'discrete_adjoint.py',
'direct_differentiation.py',
'fsi_computation.py',
- 'SU2_CFD.py'],
+ 'SU2_CFD.py'],
install_dir: join_paths(get_option('bindir')))
install_data(['SU2/__init__.py'], install_dir: join_paths(get_option('bindir'), 'SU2/'))
@@ -21,7 +21,7 @@ install_data(['SU2/__init__.py'], install_dir: join_paths(get_option('bindir'),
install_data(['SU2/eval/design.py',
'SU2/eval/functions.py',
'SU2/eval/gradients.py',
- 'SU2/eval/__init__.py'],
+ 'SU2/eval/__init__.py'],
install_dir: join_paths(get_option('bindir'), 'SU2/eval'))
install_data(['SU2/io/config.py',
@@ -32,17 +32,17 @@ install_data(['SU2/io/config.py',
'SU2/io/state.py',
'SU2/io/tools.py',
'SU2/io/historyMap.py',
- 'SU2/io/__init__.py'],
+ 'SU2/io/__init__.py'],
install_dir: join_paths(get_option('bindir'), 'SU2/io'))
install_data(['SU2/mesh/adapt.py',
'SU2/mesh/tools.py',
- 'SU2/mesh/__init__.py'],
+ 'SU2/mesh/__init__.py'],
install_dir: join_paths(get_option('bindir'), 'SU2/mesh'))
install_data(['SU2/opt/project.py',
'SU2/opt/scipy_tools.py',
- 'SU2/opt/__init__.py'],
+ 'SU2/opt/__init__.py'],
install_dir: join_paths(get_option('bindir'), 'SU2/opt'))
install_data(['SU2/run/adaptation.py',
@@ -70,14 +70,12 @@ install_data(['SU2/util/bunch.py',
'SU2/util/__init__.py'],
install_dir: join_paths(get_option('bindir'), 'SU2/util'))
-install_data(['FSI/__init__.py',
- 'FSI/FSIInterface.py',
- 'FSI/PitchPlungeAirfoilStructuralTester.py',
- 'FSI/io/__init__.py',
- 'FSI/io/FSI_config.py',
- 'FSI/util/__init__.py',
- 'FSI/util/switch.py'],
- install_dir: join_paths(get_option('bindir'), 'FSI'))
-
-
+install_data(['FSI_tools/__init__.py',
+ 'FSI_tools/FSIInterface.py',
+ 'FSI_tools/FSI_config.py',
+ 'FSI_tools/switch.py'],
+ install_dir: join_paths(get_option('bindir'), 'FSI_tools'))
+install_data(['SU2_Nastran/__init__.py',
+ 'SU2_Nastran/pysu2_nastran.py',],
+ install_dir: join_paths(get_option('bindir'), 'SU2_Nastran'))