Merge pull request #1199 from su2code/fix_legacy_deformation

Fix grid velocities in dual-time RANS simulations with deforming grids (aeroelastic and legacy FSI)

Author: Nicola-Fonzi

SU2_CFD/include/integration/CIntegration.hpp

@@ -125,20 +125,20 @@ class CIntegration {
   inline bool GetConvergence_FullMG(void) const { return Convergence_FullMG; }
 
   /*!
-   * \brief Save the solution, and volume at different time steps.
+   * \brief Save the geometry at different time steps.
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solution - Flow solution.
+   * \param[in] solver - Mesh solver.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetDualTime_Solver(CGeometry *geometry, CSolver *solver, CConfig *config, unsigned short iMesh);
+  void SetDualTime_Geometry(CGeometry *geometry, CSolver *mesh_solver, const CConfig *config, unsigned short iMesh);
 
   /*!
-   * \brief Save the structural solution at different time steps.
+   * \brief Save the solution at different time steps, and reset certain fields for the next timestep.
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver_container - Structural solution.
+   * \param[in] solver - Some solver.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetStructural_Solver(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh);
+  virtual void SetDualTime_Solver(const CGeometry *geometry, CSolver *solver, const CConfig *config, unsigned short iMesh);
 
   /*!
    * \brief A virtual member.

SU2_CFD/include/integration/CStructuralIntegration.hpp

@@ -51,6 +51,14 @@ class CStructuralIntegration final : public CIntegration {
                             CNumerics ******numerics_container, CConfig **config,
                             unsigned short RunTime_EqSystem, unsigned short iZone, unsigned short iInst) override;
 
+  /*!
+   * \brief Save the solution at different time steps, and reset certain fields for the next timestep.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver - Structural solver.
+   * \param[in] config - Definition of the problem.
+   */
+  void SetDualTime_Solver(const CGeometry *geometry, CSolver *solver, const CConfig *config, unsigned short iMesh) override;
+
 private:
   /*!
    * \brief Do the space integration of the numerical system on a FEM framework.

SU2_CFD/include/iteration/CFluidIteration.hpp

@@ -115,6 +115,8 @@ class CFluidIteration : public CIteration {
                    CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone,
                    unsigned short val_iInst) override;
 
+ private:
+
   /*!
    * \brief Imposes a gust via the grid velocities.
    * \author S. Padron
@@ -146,4 +148,11 @@ class CFluidIteration : public CIteration {
    * \return Boolean indicating weather calculation should be stopped
    */
   bool MonitorFixed_CL(COutput* output, CGeometry* geometry, CSolver** solver, CConfig* config);
+
+  /*!
+   * \brief Store old aeroelastic solutions
+   * \param[in,out] config - Definition of the particular problem.
+   */
+  void SetDualTime_Aeroelastic(CConfig* config) const;
+
 };

SU2_CFD/include/solvers/CFEASolver.hpp

@@ -314,7 +314,7 @@ class CFEASolver : public CSolver {
    * \param[in] numerics - Description of the numerical method.
    * \param[in] config - Definition of the particular problem.
    */
-  void Compute_MassMatrix(CGeometry *geometry,
+  void Compute_MassMatrix(const CGeometry *geometry,
                           CNumerics **numerics,
                           const CConfig *config) final;
 
@@ -324,7 +324,7 @@ class CFEASolver : public CSolver {
    * \param[in] numerics - Description of the numerical method.
    * \param[in] config - Definition of the particular problem.
    */
-  void Compute_MassRes(CGeometry *geometry,
+  void Compute_MassRes(const CGeometry *geometry,
                        CNumerics **numerics,
                        const CConfig *config) final;
 
@@ -472,7 +472,7 @@ class CFEASolver : public CSolver {
    * \param[in] numerics - Numerical methods.
    * \param[in] config - Definition of the particular problem.
    */
-  void ImplicitNewmark_Iteration(CGeometry *geometry,
+  void ImplicitNewmark_Iteration(const CGeometry *geometry,
                                  CNumerics **numerics,
                                  const CConfig *config) final;
 
@@ -481,22 +481,22 @@ class CFEASolver : public CSolver {
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  void ImplicitNewmark_Update(CGeometry *geometry, CConfig *config) final;
+  void ImplicitNewmark_Update(const CGeometry *geometry, const CConfig *config) final;
 
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  void ImplicitNewmark_Relaxation(CGeometry *geometry, CConfig *config) final;
+  void ImplicitNewmark_Relaxation(const CGeometry *geometry, const CConfig *config) final;
 
   /*!
    * \brief Iterate using an implicit Generalized Alpha solver.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] numerics - Numerical methods.
    * \param[in] config - Definition of the particular problem.
    */
-  void GeneralizedAlpha_Iteration(CGeometry *geometry,
+  void GeneralizedAlpha_Iteration(const CGeometry *geometry,
                                   CNumerics **numerics,
                                   const CConfig *config) final;
 
@@ -505,21 +505,21 @@ class CFEASolver : public CSolver {
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  void GeneralizedAlpha_UpdateDisp(CGeometry *geometry, CConfig *config) final;
+  void GeneralizedAlpha_UpdateDisp(const CGeometry *geometry, const CConfig *config) final;
 
   /*!
    * \brief Update the solution using an implicit Generalized Alpha solver.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  void GeneralizedAlpha_UpdateSolution(CGeometry *geometry, CConfig *config) final;
+  void GeneralizedAlpha_UpdateSolution(const CGeometry *geometry, const CConfig *config) final;
 
   /*!
    * \brief Update the solution using an implicit Generalized Alpha solver.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  void GeneralizedAlpha_UpdateLoads(CGeometry *geometry, const CConfig *config) final;
+  void GeneralizedAlpha_UpdateLoads(const CGeometry *geometry, const CConfig *config) final;
 
   /*!
    * \brief Postprocessing.

SU2_CFD/include/solvers/CSolver.hpp

@@ -1557,7 +1557,7 @@ class CSolver {
    * \param[in] numerics - Numerical methods.
    * \param[in] config - Definition of the particular problem.
    */
-  inline virtual void ImplicitNewmark_Iteration(CGeometry *geometry,
+  inline virtual void ImplicitNewmark_Iteration(const CGeometry *geometry,
                                                 CNumerics **numerics,
                                                 const CConfig *config) { }
 
@@ -1567,25 +1567,25 @@ class CSolver {
    * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] config - Definition of the particular problem.
    */
-  inline virtual void ImplicitNewmark_Update(CGeometry *geometry,
-                                             CConfig *config) { }
+  inline virtual void ImplicitNewmark_Update(const CGeometry *geometry,
+                                             const CConfig *config) { }
 
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] config - Definition of the particular problem.
    */
-  inline virtual void ImplicitNewmark_Relaxation(CGeometry *geometry,
-                                                 CConfig *config) { }
+  inline virtual void ImplicitNewmark_Relaxation(const CGeometry *geometry,
+                                                 const CConfig *config) { }
 
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] numerics - Numerical methods.
    * \param[in] config - Definition of the particular problem.
    */
-  inline virtual void GeneralizedAlpha_Iteration(CGeometry *geometry,
+  inline virtual void GeneralizedAlpha_Iteration(const CGeometry *geometry,
                                                  CNumerics **numerics,
                                                  const CConfig *config) { }
 
@@ -1595,25 +1595,25 @@ class CSolver {
    * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] config - Definition of the particular problem.
    */
-  inline virtual void GeneralizedAlpha_UpdateDisp(CGeometry *geometry,
-                                                  CConfig *config) { }
+  inline virtual void GeneralizedAlpha_UpdateDisp(const CGeometry *geometry,
+                                                  const CConfig *config) { }
 
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] config - Definition of the particular problem.
    */
-  inline virtual void GeneralizedAlpha_UpdateSolution(CGeometry *geometry,
-                                                      CConfig *config) { }
+  inline virtual void GeneralizedAlpha_UpdateSolution(const CGeometry *geometry,
+                                                      const CConfig *config) { }
 
   /*!
    * \brief A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] solver_container - Container vector with all the solutions.
    * \param[in] config - Definition of the particular problem.
    */
-  inline virtual void GeneralizedAlpha_UpdateLoads(CGeometry *geometry,
+  inline virtual void GeneralizedAlpha_UpdateLoads(const CGeometry *geometry,
                                                    const CConfig *config) { }
 
   /*!
@@ -3770,7 +3770,7 @@ class CSolver {
    * \param[in] numerics - Description of the numerical method.
    * \param[in] config - Definition of the particular problem.
    */
-  inline virtual void Compute_MassMatrix(CGeometry *geometry,
+  inline virtual void Compute_MassMatrix(const CGeometry *geometry,
                                          CNumerics **numerics,
                                          const CConfig *config) { }
 
@@ -3780,7 +3780,7 @@ class CSolver {
    * \param[in] numerics - Description of the numerical method.
    * \param[in] config - Definition of the particular problem.
    */
-  inline virtual void Compute_MassRes(CGeometry *geometry,
+  inline virtual void Compute_MassRes(const CGeometry *geometry,
                                       CNumerics **numerics,
                                       const CConfig *config) { }
 

SU2_CFD/src/integration/CIntegration.cpp

@@ -208,14 +208,10 @@ void CIntegration::Time_Integration(CGeometry *geometry, CSolver **solver_contai
 
 }
 
-void CIntegration::SetDualTime_Solver(CGeometry *geometry, CSolver *solver, CConfig *config, unsigned short iMesh) {
+void CIntegration::SetDualTime_Geometry(CGeometry *geometry, CSolver *mesh_solver, const CConfig *config, unsigned short iMesh) {
 
   SU2_OMP_PARALLEL
   {
-  /*--- Store old solution, volumes and coordinates (in case there is grid movement). ---*/
-
-  solver->GetNodes()->Set_Solution_time_n1();
-  solver->GetNodes()->Set_Solution_time_n();
 
   geometry->nodes->SetVolume_nM1();
   geometry->nodes->SetVolume_n();
@@ -225,6 +221,19 @@ void CIntegration::SetDualTime_Solver(CGeometry *geometry, CSolver *solver, CCon
     geometry->nodes->SetCoord_n();
   }
 
+  if ((iMesh==MESH_0) && config->GetDeform_Mesh()) mesh_solver->SetDualTime_Mesh();
+
+  } // end SU2_OMP_PARALLEL
+}
+
+void CIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolver *solver, const CConfig *config, unsigned short iMesh) {
+
+  SU2_OMP_PARALLEL
+  {
+  /*--- Store old solution ---*/
+  solver->GetNodes()->Set_Solution_time_n1();
+  solver->GetNodes()->Set_Solution_time_n();
+
   SU2_OMP_MASTER
   solver->ResetCFLAdapt();
   SU2_OMP_BARRIER
@@ -239,113 +248,5 @@ void CIntegration::SetDualTime_Solver(CGeometry *geometry, CSolver *solver, CCon
     solver->GetNodes()->SetLocalCFL(iPoint, config->GetCFL(iMesh));
   }
 
-  /*--- Store old aeroelastic solutions ---*/
-  SU2_OMP_MASTER
-  if (config->GetGrid_Movement() && config->GetAeroelastic_Simulation() && (iMesh == MESH_0)) {
-
-    config->SetAeroelastic_n1();
-    config->SetAeroelastic_n();
-
-    /*--- Also communicate plunge and pitch to the master node. Needed for output in case of parallel run ---*/
-
-#ifdef HAVE_MPI
-    su2double plunge, pitch, *plunge_all = NULL, *pitch_all = NULL;
-    unsigned short iMarker, iMarker_Monitoring;
-    unsigned long iProcessor, owner, *owner_all = NULL;
-
-    string Marker_Tag, Monitoring_Tag;
-    int nProcessor = size;
-
-    /*--- Only if master node allocate memory ---*/
-
-    if (rank == MASTER_NODE) {
-      plunge_all = new su2double[nProcessor];
-      pitch_all  = new su2double[nProcessor];
-      owner_all  = new unsigned long[nProcessor];
-    }
-
-    /*--- Find marker and give it's plunge and pitch coordinate to the master node ---*/
-
-    for (iMarker_Monitoring = 0; iMarker_Monitoring < config->GetnMarker_Monitoring(); iMarker_Monitoring++) {
-
-      for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
-
-        Monitoring_Tag = config->GetMarker_Monitoring_TagBound(iMarker_Monitoring);
-        Marker_Tag = config->GetMarker_All_TagBound(iMarker);
-        if (Marker_Tag == Monitoring_Tag) { owner = 1; break;
-        } else {
-          owner = 0;
-        }
-
-      }
-      plunge = config->GetAeroelastic_plunge(iMarker_Monitoring);
-      pitch  = config->GetAeroelastic_pitch(iMarker_Monitoring);
-
-      /*--- Gather the data on the master node. ---*/
-
-      SU2_MPI::Gather(&plunge, 1, MPI_DOUBLE, plunge_all, 1, MPI_DOUBLE, MASTER_NODE, SU2_MPI::GetComm());
-      SU2_MPI::Gather(&pitch, 1, MPI_DOUBLE, pitch_all, 1, MPI_DOUBLE, MASTER_NODE, SU2_MPI::GetComm());
-      SU2_MPI::Gather(&owner, 1, MPI_UNSIGNED_LONG, owner_all, 1, MPI_UNSIGNED_LONG, MASTER_NODE, SU2_MPI::GetComm());
-
-      /*--- Set plunge and pitch on the master node ---*/
-
-      if (rank == MASTER_NODE) {
-        for (iProcessor = 0; iProcessor < (unsigned long)nProcessor; iProcessor++) {
-          if (owner_all[iProcessor] == 1) {
-            config->SetAeroelastic_plunge(iMarker_Monitoring, plunge_all[iProcessor]);
-            config->SetAeroelastic_pitch(iMarker_Monitoring, pitch_all[iProcessor]);
-            break;
-          }
-        }
-      }
-
-    }
-
-    if (rank == MASTER_NODE) {
-      delete [] plunge_all;
-      delete [] pitch_all;
-      delete [] owner_all;
-    }
-#endif
-  }
-  SU2_OMP_BARRIER
-
   } // end SU2_OMP_PARALLEL
-
-}
-
-void CIntegration::SetStructural_Solver(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh) {
-
-  bool fsi = config->GetFSI_Simulation();
-
-  /*--- Update the solution according to the integration scheme used ---*/
-
-  switch (config->GetKind_TimeIntScheme_FEA()) {
-    case (CD_EXPLICIT):
-      break;
-    case (NEWMARK_IMPLICIT):
-      if (fsi) solver_container[FEA_SOL]->ImplicitNewmark_Relaxation(geometry, config);
-      break;
-    case (GENERALIZED_ALPHA):
-      solver_container[FEA_SOL]->GeneralizedAlpha_UpdateSolution(geometry, config);
-      solver_container[FEA_SOL]->GeneralizedAlpha_UpdateLoads(geometry, config);
-      break;
-  }
-
-  /*--- Store the solution at t+1 as solution at t, both for the local points and for the halo points ---*/
-
-  solver_container[FEA_SOL]->GetNodes()->Set_Solution_time_n();
-  solver_container[FEA_SOL]->GetNodes()->SetSolution_Vel_time_n();
-  solver_container[FEA_SOL]->GetNodes()->SetSolution_Accel_time_n();
-
-  /*--- If FSI problem, save the last Aitken relaxation parameter of the previous time step ---*/
-
-  if (fsi) {
-
-    su2double WAitk=0.0;
-
-    WAitk = solver_container[FEA_SOL]->GetWAitken_Dyn();
-    solver_container[FEA_SOL]->SetWAitken_Dyn_tn1(WAitk);
-
-  }
 }