diff --git a/SU2_CFD/include/integration/CIntegration.hpp b/SU2_CFD/include/integration/CIntegration.hpp
index 0c88a5d05596..2062713e718a 100644
--- a/SU2_CFD/include/integration/CIntegration.hpp
+++ b/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.
diff --git a/SU2_CFD/include/integration/CStructuralIntegration.hpp b/SU2_CFD/include/integration/CStructuralIntegration.hpp
index 845167a3d4b2..bcd446935803 100644
--- a/SU2_CFD/include/integration/CStructuralIntegration.hpp
+++ b/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.
diff --git a/SU2_CFD/include/iteration/CFluidIteration.hpp b/SU2_CFD/include/iteration/CFluidIteration.hpp
index 82cd95628be0..9271ef475a16 100644
--- a/SU2_CFD/include/iteration/CFluidIteration.hpp
+++ b/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;
+
 };
diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp
index d845127be60d..2e3cee9e5e09 100644
--- a/SU2_CFD/include/solvers/CFEASolver.hpp
+++ b/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,14 +481,14 @@ 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.
@@ -496,7 +496,7 @@ class CFEASolver : public CSolver {
    * \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.
diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp
index 062711cb47b6..5b11ab8b4208 100644
--- a/SU2_CFD/include/solvers/CSolver.hpp
+++ b/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,8 +1567,8 @@ 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.
@@ -1576,8 +1576,8 @@ class CSolver {
    * \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.
@@ -1585,7 +1585,7 @@ class CSolver {
    * \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,8 +1595,8 @@ 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.
@@ -1604,8 +1604,8 @@ class CSolver {
    * \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.
@@ -1613,7 +1613,7 @@ class CSolver {
    * \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) { }
 
diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp
index 811b1b5608cd..c2699c9eb560 100644
--- a/SU2_CFD/src/integration/CIntegration.cpp
+++ b/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);
-
-  }
 }
diff --git a/SU2_CFD/src/integration/CStructuralIntegration.cpp b/SU2_CFD/src/integration/CStructuralIntegration.cpp
index 3667eba307a4..52fb4177ceb8 100644
--- a/SU2_CFD/src/integration/CStructuralIntegration.cpp
+++ b/SU2_CFD/src/integration/CStructuralIntegration.cpp
@@ -190,3 +190,39 @@ void CStructuralIntegration::Time_Integration_FEM(CGeometry *geometry, CSolver *
   }
 
 }
+
+void CStructuralIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolver *solver, const 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->ImplicitNewmark_Relaxation(geometry, config);
+      break;
+    case (GENERALIZED_ALPHA):
+      solver->GeneralizedAlpha_UpdateSolution(geometry, config);
+      solver->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->GetNodes()->Set_Solution_time_n();
+  solver->GetNodes()->SetSolution_Vel_time_n();
+  solver->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->GetWAitken_Dyn();
+    solver->SetWAitken_Dyn_tn1(WAitk);
+
+  }
+}
diff --git a/SU2_CFD/src/iteration/CAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CAdjFluidIteration.cpp
index 315b0d40be6c..870d6b00d2d8 100644
--- a/SU2_CFD/src/iteration/CAdjFluidIteration.cpp
+++ b/SU2_CFD/src/iteration/CAdjFluidIteration.cpp
@@ -181,6 +181,11 @@ void CAdjFluidIteration::Update(COutput* output, CIntegration**** integration, C
       integration[val_iZone][val_iInst][ADJFLOW_SOL]->SetDualTime_Solver(
           geometry[val_iZone][val_iInst][iMesh], solver[val_iZone][val_iInst][iMesh][ADJFLOW_SOL], config[val_iZone],
           iMesh);
+
+      integration[val_iZone][val_iInst][ADJFLOW_SOL]->SetDualTime_Geometry(
+          geometry[val_iZone][val_iInst][iMesh], solver[val_iZone][val_iInst][iMesh][MESH_SOL], config[val_iZone],
+          iMesh);
+
       integration[val_iZone][val_iInst][ADJFLOW_SOL]->SetConvergence(false);
     }
 
diff --git a/SU2_CFD/src/iteration/CFEAIteration.cpp b/SU2_CFD/src/iteration/CFEAIteration.cpp
index 11b553134d61..1263c41b36e2 100644
--- a/SU2_CFD/src/iteration/CFEAIteration.cpp
+++ b/SU2_CFD/src/iteration/CFEAIteration.cpp
@@ -196,8 +196,9 @@ void CFEAIteration::Update(COutput* output, CIntegration**** integration, CGeome
   /*----------------- Update structural solver ----------------------*/
 
   if (dynamic) {
-    integration[val_iZone][val_iInst][FEA_SOL]->SetStructural_Solver(
-        geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0], config[val_iZone], MESH_0);
+    integration[val_iZone][val_iInst][FEA_SOL]->SetDualTime_Solver(
+        geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0][FEA_SOL], config[val_iZone],
+        MESH_0);
     integration[val_iZone][val_iInst][FEA_SOL]->SetConvergence(false);
 
     /*--- Verify convergence criteria (based on total time) ---*/
diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp
index 47a6f2625bb9..539252c76712 100644
--- a/SU2_CFD/src/iteration/CFluidIteration.cpp
+++ b/SU2_CFD/src/iteration/CFluidIteration.cpp
@@ -172,13 +172,15 @@ void CFluidIteration::Update(COutput* output, CIntegration**** integration, CGeo
       integration[val_iZone][val_iInst][FLOW_SOL]->SetDualTime_Solver(geometry[val_iZone][val_iInst][iMesh],
                                                                       solver[val_iZone][val_iInst][iMesh][FLOW_SOL],
                                                                       config[val_iZone], iMesh);
+
+      integration[val_iZone][val_iInst][FLOW_SOL]->SetDualTime_Geometry(geometry[val_iZone][val_iInst][iMesh],
+                                                                        solver[val_iZone][val_iInst][iMesh][MESH_SOL],
+                                                                        config[val_iZone], iMesh);
+
       integration[val_iZone][val_iInst][FLOW_SOL]->SetConvergence(false);
     }
 
-    /*--- Update dual time solver for the dynamic mesh solver ---*/
-    if (config[val_iZone]->GetDeform_Mesh()) {
-      solver[val_iZone][val_iInst][MESH_0][MESH_SOL]->SetDualTime_Mesh();
-    }
+    SetDualTime_Aeroelastic(config[val_iZone]);
 
     /*--- Update dual time solver for the turbulence model ---*/
 
@@ -199,6 +201,15 @@ void CFluidIteration::Update(COutput* output, CIntegration**** integration, CGeo
                                                                        config[val_iZone], MESH_0);
       integration[val_iZone][val_iInst][TRANS_SOL]->SetConvergence(false);
     }
+
+    /*--- Update dual time solver for the weakly coupled energy equation ---*/
+
+    if (config[val_iZone]->GetWeakly_Coupled_Heat()) {
+      integration[val_iZone][val_iInst][HEAT_SOL]->SetDualTime_Solver(geometry[val_iZone][val_iInst][MESH_0],
+                                                                      solver[val_iZone][val_iInst][MESH_0][HEAT_SOL],
+                                                                      config[val_iZone], MESH_0);
+      integration[val_iZone][val_iInst][HEAT_SOL]->SetConvergence(false);
+    }
   }
 }
 
@@ -557,3 +568,74 @@ bool CFluidIteration::MonitorFixed_CL(COutput *output, CGeometry *geometry, CSol
   /* --- Set convergence based on fixed CL convergence  --- */
   return fixed_cl_convergence;
 }
+
+void CFluidIteration::SetDualTime_Aeroelastic(CConfig* config) const {
+
+  /*--- Store old aeroelastic solutions ---*/
+
+  if (config->GetGrid_Movement() && config->GetAeroelastic_Simulation()) {
+
+    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 = nullptr, *pitch_all = nullptr;
+    unsigned short iMarker, iMarker_Monitoring;
+    unsigned long iProcessor, owner, *owner_all = nullptr;
+
+    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;
+          }
+        }
+      }
+    }
+
+    delete [] plunge_all;
+    delete [] pitch_all;
+    delete [] owner_all;
+#endif
+  }
+
+}
diff --git a/SU2_CFD/src/iteration/CHeatIteration.cpp b/SU2_CFD/src/iteration/CHeatIteration.cpp
index dbdd69ce22d3..adf02634aa18 100644
--- a/SU2_CFD/src/iteration/CHeatIteration.cpp
+++ b/SU2_CFD/src/iteration/CHeatIteration.cpp
@@ -56,6 +56,11 @@ void CHeatIteration::Update(COutput* output, CIntegration**** integration, CGeom
       integration[val_iZone][val_iInst][HEAT_SOL]->SetDualTime_Solver(geometry[val_iZone][val_iInst][iMesh],
                                                                       solver[val_iZone][val_iInst][iMesh][HEAT_SOL],
                                                                       config[val_iZone], iMesh);
+
+      integration[val_iZone][val_iInst][HEAT_SOL]->SetDualTime_Geometry(geometry[val_iZone][val_iInst][iMesh],
+                                                                        solver[val_iZone][val_iInst][iMesh][MESH_SOL],
+                                                                        config[val_iZone], iMesh);
+
       integration[val_iZone][val_iInst][HEAT_SOL]->SetConvergence(false);
     }
   }
diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp
index c54b3bb9706d..66fe4d41ba58 100644
--- a/SU2_CFD/src/solvers/CFEASolver.cpp
+++ b/SU2_CFD/src/solvers/CFEASolver.cpp
@@ -936,7 +936,7 @@ void CFEASolver::Compute_StiffMatrix_NodalStressRes(CGeometry *geometry, CNumeri
 
 }
 
-void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
+void CFEASolver::Compute_MassMatrix(const CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
 
   const bool topology_mode = config->GetTopology_Optimization();
   const su2double simp_minstiff = config->GetSIMP_MinStiffness();
@@ -1021,7 +1021,7 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, c
 
 }
 
-void CFEASolver::Compute_MassRes(CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
+void CFEASolver::Compute_MassRes(const CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
 
   const bool topology_mode = config->GetTopology_Optimization();
   const su2double simp_minstiff = config->GetSIMP_MinStiffness();
@@ -2183,7 +2183,7 @@ su2double CFEASolver::Compute_LoadCoefficient(su2double CurrentTime, su2double R
 
 }
 
-void CFEASolver::ImplicitNewmark_Iteration(CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
+void CFEASolver::ImplicitNewmark_Iteration(const CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
 
   const bool first_iter = (config->GetInnerIter() == 0);
   const bool dynamic = (config->GetTime_Domain());
@@ -2267,7 +2267,7 @@ void CFEASolver::ImplicitNewmark_Iteration(CGeometry *geometry, CNumerics **nume
 
 }
 
-void CFEASolver::ImplicitNewmark_Update(CGeometry *geometry, CConfig *config) {
+void CFEASolver::ImplicitNewmark_Update(const CGeometry *geometry, const CConfig *config) {
 
   const bool dynamic = (config->GetTime_Domain());
 
@@ -2314,7 +2314,7 @@ void CFEASolver::ImplicitNewmark_Update(CGeometry *geometry, CConfig *config) {
   } // end SU2_OMP_PARALLEL
 }
 
-void CFEASolver::ImplicitNewmark_Relaxation(CGeometry *geometry, CConfig *config) {
+void CFEASolver::ImplicitNewmark_Relaxation(const CGeometry *geometry, const CConfig *config) {
 
   const bool dynamic = (config->GetTime_Domain());
 
@@ -2362,7 +2362,7 @@ void CFEASolver::ImplicitNewmark_Relaxation(CGeometry *geometry, CConfig *config
 }
 
 
-void CFEASolver::GeneralizedAlpha_Iteration(CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
+void CFEASolver::GeneralizedAlpha_Iteration(const CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
 
   const bool first_iter = (config->GetInnerIter() == 0);
   const bool dynamic = (config->GetTime_Domain());
@@ -2470,7 +2470,7 @@ void CFEASolver::GeneralizedAlpha_Iteration(CGeometry *geometry, CNumerics **num
 
 }
 
-void CFEASolver::GeneralizedAlpha_UpdateDisp(CGeometry *geometry, CConfig *config) {
+void CFEASolver::GeneralizedAlpha_UpdateDisp(const CGeometry *geometry, const CConfig *config) {
 
   /*--- Update displacement components of the solution. ---*/
 
@@ -2481,7 +2481,7 @@ void CFEASolver::GeneralizedAlpha_UpdateDisp(CGeometry *geometry, CConfig *confi
 
 }
 
-void CFEASolver::GeneralizedAlpha_UpdateSolution(CGeometry *geometry, CConfig *config) {
+void CFEASolver::GeneralizedAlpha_UpdateSolution(const CGeometry *geometry, const CConfig *config) {
 
   const su2double alpha_f = config->Get_Int_Coeffs(2);
   const su2double alpha_m = config->Get_Int_Coeffs(3);
@@ -2535,7 +2535,7 @@ void CFEASolver::GeneralizedAlpha_UpdateSolution(CGeometry *geometry, CConfig *c
 
 }
 
-void CFEASolver::GeneralizedAlpha_UpdateLoads(CGeometry *geometry, const CConfig *config) {
+void CFEASolver::GeneralizedAlpha_UpdateLoads(const CGeometry *geometry, const CConfig *config) {
 
   /*--- Set the load conditions of the time step n+1 as the load conditions for time step n ---*/
   nodes->Set_SurfaceLoad_Res_n();
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index e14d554ca9e0..6a9fbe136f14 100644
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -1313,7 +1313,7 @@ def main():
     pywrapper_rigidMotion.cfg_dir       = "py_wrapper/flatPlate_rigidMotion"
     pywrapper_rigidMotion.cfg_file      = "flatPlate_rigidMotion_Conf.cfg"
     pywrapper_rigidMotion.test_iter     = 5
-    pywrapper_rigidMotion.test_vals     = [-1.614165, 2.242641, -0.038307, 0.173866]
+    pywrapper_rigidMotion.test_vals     = [-1.614164, 2.242568, -0.028488, 0.173947]
     pywrapper_rigidMotion.su2_exec      = "mpirun -np 2 python launch_flatPlate_rigidMotion.py --parallel -f"
     pywrapper_rigidMotion.timeout       = 1600
     pywrapper_rigidMotion.tol           = 0.00001
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index b99d15288ebf..bd893681abb0 100644
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -1892,7 +1892,7 @@ def main():
     pywrapper_rigidMotion.cfg_dir       = "py_wrapper/flatPlate_rigidMotion"
     pywrapper_rigidMotion.cfg_file      = "flatPlate_rigidMotion_Conf.cfg"
     pywrapper_rigidMotion.test_iter     = 5
-    pywrapper_rigidMotion.test_vals     = [-1.614167, 2.242632, -0.037871, 0.173912]
+    pywrapper_rigidMotion.test_vals     = [-1.614167, 2.242558, -0.027574, 0.173990]
     pywrapper_rigidMotion.su2_exec      = "python launch_flatPlate_rigidMotion.py -f"
     pywrapper_rigidMotion.new_output      = True
     pywrapper_rigidMotion.timeout       = 1600