[WIP] Change for sliding plane when relative frame is used

Common/src/interface_interpolation/CNearestNeighbor.cpp

@@ -117,7 +117,55 @@ void CNearestNeighbor::SetTransferCoeff(const CConfig* const* config) {
             const auto idx = iProcessor * MaxLocalVertex_Donor + jVertex;
             const auto pGlobalPoint = Buffer_Receive_GlobalPoint[idx];
             const su2double* Coord_j = Buffer_Receive_Coord[idx];
-            const auto dist2 = GeometryToolbox::SquaredDistance(nDim, Coord_i, Coord_j);
+
+            su2double Theta, Phi, Psi;
+            su2double rotCoord_i[3] = {0.0, 0.0, 0.0}, rotCoord_j[3] = {0.0, 0.0, 0.0};
+            su2double rotMatrix[3][3] = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
+            const su2double zeros[3] = {0.0};
+            for (unsigned short iDim = 0; iDim < 3; iDim++) {
+              rotCoord_i[iDim] = Coord_i[iDim];
+              rotCoord_j[iDim] = Coord_j[iDim];
+            }
+
+            if (config[targetZone]->GetRotating_Frame() == YES) {
+              su2double Omega_i[3] = {0.0, 0.0, 0.0};
+              su2double dt = config[targetZone]->GetDelta_UnstTimeND();
+              unsigned long TimeIter = config[targetZone]->GetTimeIter();
+              for (unsigned short iDim = 0; iDim < 3; iDim++) {
+                Omega_i[iDim] = config[targetZone]->GetRotation_Rate(iDim) / config[targetZone]->GetOmega_Ref();
+              }
+
+              /*--- Compute the rotation matrix. Note that the implicit
+              ordering is rotation about the x-axis, y-axis, then z-axis. ---*/
+              Theta = Omega_i[0] * dt * TimeIter;
+              Phi = Omega_i[1] * dt * TimeIter;
+              Psi = Omega_i[2] * dt * TimeIter;
+              GeometryToolbox::RotationMatrix(Theta, Phi, Psi, rotMatrix);
+
+              /*--- Compute transformed point coordinates. ---*/
+              GeometryToolbox::Rotate(rotMatrix, zeros, Coord_i, rotCoord_i);
+            }
+
+            if (config[donorZone]->GetRotating_Frame() == YES) {
+              su2double Omega_j[3] = {0.0, 0.0, 0.0};
+              su2double dt = config[donorZone]->GetDelta_UnstTimeND();
+              unsigned long TimeIter = config[donorZone]->GetTimeIter();
+              for (unsigned short iDim = 0; iDim < 3; iDim++) {
+                Omega_j[iDim] = config[donorZone]->GetRotation_Rate(iDim) / config[donorZone]->GetOmega_Ref();
+              }
+
+              /*--- Compute the rotation matrix. Note that the implicit
+              ordering is rotation about the x-axis, y-axis, then z-axis. ---*/
+              Theta = Omega_j[0] * dt * TimeIter;
+              Phi = Omega_j[1] * dt * TimeIter;
+              Psi = Omega_j[2] * dt * TimeIter;
+              GeometryToolbox::RotationMatrix(Theta, Phi, Psi, rotMatrix);
+
+              /*--- Compute transformed point coordinates. ---*/
+              GeometryToolbox::Rotate(rotMatrix, zeros, Coord_j, rotCoord_j);
+            }
+
+            const auto dist2 = GeometryToolbox::SquaredDistance(nDim, rotCoord_i, rotCoord_j);
 
             donorInfo[iDonor++] = DonorInfo(dist2, pGlobalPoint, iProcessor);
           }

SU2_CFD/include/interfaces/CInterface.hpp

@@ -125,6 +125,16 @@ class CInterface {
                                  const CConfig *donor_config, unsigned long Marker_Donor,
                                  unsigned long Vertex_Donor, unsigned long Point_Donor) = 0;
 
+  /*!
+   * \brief A virtual member.
+   * \param[in] donor_config - Definition of the problem at the donor mesh.   
+   * \param[in] donor_geometry - Geometry of the donor mesh.
+   * \param[in] target_config - Definition of the problem at the donor mesh.   
+   * \param[in] target_geometry - Geometry of the donor mesh.
+   */
+  inline virtual void GetDonor_Velocity_RotatingFrame(const CConfig *donor_config, CGeometry *donor_geometry, 
+                                 const CConfig *target_config, CGeometry *target_geometry) {}
+
   /*!
    * \brief Initializes the target variable.
    * \param[in] target_solution - Solution from the target mesh.

SU2_CFD/include/interfaces/cfd/CSlidingInterface.hpp

@@ -55,6 +55,17 @@ class CSlidingInterface : public CInterface {
                          unsigned long Marker_Donor, unsigned long Vertex_Donor, unsigned long Point_Donor) override;
 
   /*!
+  * \brief Rotate the velocity if rotating frame is applied.
+   * \param[in] donor_config - Definition of the problem at the donor mesh.   
+   * \param[in] donor_geometry - Geometry of the donor mesh.
+   * \param[in] target_config - Definition of the problem at the donor mesh.   
+   * \param[in] target_geometry - Geometry of the donor mesh.
+   */
+  void GetDonor_Velocity_RotatingFrame(const CConfig *donor_config, CGeometry *donor_geometry, 
+                                 const CConfig *target_config, CGeometry *target_geometry) override;
+
+
+  /*!  
    * \brief A virtual member, initializes the target variable for sliding mesh.
    * \param[in] target_solution - Solution from the target mesh.
    * \param[in] Marker_Target - Index of the target marker.

SU2_CFD/src/drivers/CMultizoneDriver.cpp

@@ -101,7 +101,7 @@ CMultizoneDriver::CMultizoneDriver(char* confFile, unsigned short val_nZone, SU2
   prefixed_motion = new bool[nZone];
   for (iZone = 0; iZone < nZone; iZone++){
     switch (config_container[iZone]->GetKind_GridMovement()){
-      case RIGID_MOTION:
+      case RIGID_MOTION: case ROTATING_FRAME:
         prefixed_motion[iZone] = true; break;
       default:
         prefixed_motion[iZone] = false; break;
@@ -269,7 +269,7 @@ void CMultizoneDriver::Preprocess(unsigned long TimeIter) {
   if (driver_config->GetTime_Domain()) {
     for (iZone = 0; iZone < nZone; iZone++) {
       for (unsigned short jZone = 0; jZone < nZone; jZone++){
-        if(jZone != iZone && interpolator_container[iZone][jZone] != nullptr && prefixed_motion[iZone])
+        if(jZone != iZone && interpolator_container[iZone][jZone] != nullptr && (prefixed_motion[iZone] || prefixed_motion[jZone]))
           interpolator_container[iZone][jZone]->SetTransferCoeff(config_container);
       }
     }

SU2_CFD/src/interfaces/CInterface.cpp

@@ -115,6 +115,12 @@ void CInterface::BroadcastData(const CInterpolator& interpolator,
         if (donor_geometry->nodes->GetDomain(iPoint)) {
 
           GetDonor_Variable(donor_solution, donor_geometry, donor_config, markDonor, iVertex, iPoint);
+
+          /*---rotate the velocity if rotating frame is applied---*/
+          if (donor_solution->GetnPrimVar() > 2){
+            GetDonor_Velocity_RotatingFrame(donor_config, donor_geometry, target_config, target_geometry);
+          }
+
           for (auto iVar = 0u; iVar < nVar; iVar++) sendDonorVar(iSend, iVar) = Donor_Variable[iVar];
 
           sendDonorIdx[iSend] = donor_geometry->nodes->GetGlobalIndex(iPoint);

SU2_CFD/src/interfaces/cfd/CSlidingInterface.cpp

@@ -30,6 +30,7 @@
 #include "../../../../Common/include/CConfig.hpp"
 #include "../../../../Common/include/geometry/CGeometry.hpp"
 #include "../../../include/solvers/CSolver.hpp"
+#include "../../../../Common/include/toolboxes/geometry_toolbox.hpp"
 
 CSlidingInterface::CSlidingInterface(unsigned short val_nVar, unsigned short val_nConst) : CInterface() {
 
@@ -43,6 +44,81 @@
 
 }
 
+void CSlidingInterface::GetDonor_Velocity_RotatingFrame(const CConfig *donor_config, CGeometry *donor_geometry, 
+                                 const CConfig *target_config, CGeometry *target_geometry){
+
+  /*---rotate the velocity for rotating frame---*/
+  if (donor_config->GetRotating_Frame()==YES){
+
+    unsigned short nDim = donor_geometry->GetnDim();
+
+    su2double Theta, Phi, Psi;
+    su2double oriVel[3] = {0.0, 0.0, 0.0};
+    su2double rotVel[3] = {0.0, 0.0, 0.0};
+    su2double rotMatrix[3][3] = {{1.0,0.0,0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}};
+    su2double Omega[3] = {0.0, 0.0, 0.0};
+    const su2double zeros[3] = {0.0};
+    su2double dt = donor_config->GetDelta_UnstTimeND();
+    unsigned long TimeIter = donor_config->GetTimeIter();
+    for (unsigned short iDim=0; iDim<3; iDim++){
+      Omega[iDim] = donor_config->GetRotation_Rate(iDim)/donor_config->GetOmega_Ref();
+    }
+
+    /*--- Compute the rotation matrix. Note that the implicit
+     ordering is rotation about the x-axis, y-axis, then z-axis. ---*/
+    Theta    = Omega[0]*dt*TimeIter;      Phi = Omega[1]*dt*TimeIter;     Psi = Omega[2]*dt*TimeIter;
+    GeometryToolbox::RotationMatrix(Theta, Phi, Psi, rotMatrix);
+
+    /*--- Velocities before rotating ---*/
+    for (unsigned short iDim = 0; iDim < nDim; ++iDim)
+      oriVel[iDim] = Donor_Variable[iDim+1];
+
+    /*--- Compute transformed velocities. ---*/
+    GeometryToolbox::Rotate(rotMatrix, zeros, oriVel, rotVel);
+
+    /*--- set the rotated velocity ---*/
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      Donor_Variable[iDim+1] = rotVel[iDim];
+    }
+
+  }
+
+  if (target_config->GetRotating_Frame()==YES){
+
+    unsigned short nDim = target_geometry->GetnDim();
+
+    su2double Theta, Phi, Psi, cosTheta, sinTheta, cosPhi, sinPhi, cosPsi, sinPsi;
+    su2double oriVel[3] = {0.0, 0.0, 0.0};
+    su2double rotVel[3] = {0.0, 0.0, 0.0};
+    su2double rotMatrix[3][3] = {{1.0,0.0,0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}};
+    su2double Omega[3] = {0.0, 0.0, 0.0};
+    const su2double zeros[3] = {0.0};
+    su2double dt = target_config->GetDelta_UnstTimeND();
+    unsigned long TimeIter = donor_config->GetTimeIter();
+    for (unsigned short iDim=0; iDim<3; iDim++){
+      Omega[iDim] = -target_config->GetRotation_Rate(iDim)/target_config->GetOmega_Ref();
+    }
+
+    /*--- Compute the rotation matrix. Note that the implicit
+     ordering is rotation about the x-axis, y-axis, then z-axis. ---*/
+    Theta    = Omega[0]*dt*TimeIter;      Phi = Omega[1]*dt*TimeIter;     Psi = Omega[2]*dt*TimeIter;
+    GeometryToolbox::RotationMatrix(Theta, Phi, Psi, rotMatrix);
+
+    /*--- velocities before rotating ---*/
+    for (unsigned short iDim = 0; iDim < nDim; ++iDim)
+      oriVel[iDim] = Donor_Variable[iDim+1];
+
+    /*--- Compute transformed velocities. ---*/
+    GeometryToolbox::Rotate(rotMatrix, zeros, oriVel, rotVel);
+
+    /*--- set the rotated velocity ---*/
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      Donor_Variable[iDim+1] = rotVel[iDim];
+    }
+
+  }
+
+}
 void CSlidingInterface::GetDonor_Variable(CSolver *donor_solution, CGeometry *donor_geometry,
                                           const CConfig *donor_config, unsigned long Marker_Donor,
                                           unsigned long Vertex_Donor, unsigned long Point_Donor) {

TestCases/turbomachinery/axial_stage_full_annulus_unsteady/rotor_FA_Unsteady.cfg

@@ -0,0 +1,22 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: 2D Axial stage                                             %
+% Author: S. Vitale                                                            %
+% Institution: Delft University of Technology                                  %
+% Date: Feb 28th, 2017                                                         %
+% File Version 5.0.0"cardinal"                                                 %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% ----------------------- DYNAMIC MESH DEFINITION -----------------------------%
+%
+% Type of dynamic mesh (NONE, ROTATING_FRAME)
+GRID_MOVEMENT= ROTATING_FRAME
+%
+%
+% Angular velocity vector (rad/s) about the motion origi. Example 1250 RPM -> 130.89969389957471 rad/s
+ROTATION_RATE= 0.0 0.0 -1800.0
+%
+% Mesh input file
+MESH_FILENAME= rotor35_Full.su2