diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index 6e504444b190..f6725886df8c 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/Common/include/geometry/dual_grid/CPoint.hpp @@ -60,6 +60,9 @@ class CPoint { su2activevector Volume; /*!< \brief Volume or Area of the control volume in 3D and 2D. */ su2activevector Volume_n; /*!< \brief Volume at time n. */ su2activevector Volume_nM1; /*!< \brief Volume at time n-1. */ + su2activevector Volume_Old; /*!< \brief Old containers for Volume. */ + su2activevector Volume_n_Old; /*!< \brief Old containers for Volume at time n. */ + su2activevector Volume_nM1_Old; /*!< \brief Old containers for Volume at time n-1. */ su2activevector Periodic_Volume; /*!< \brief Missing component of volume or area of a control volume on a periodic marker in 3D and 2D. */ su2vector Domain; /*!< \brief Indicates if a point must be computed or belong to another boundary */ @@ -524,9 +527,9 @@ class CPoint { inline su2double GetVolume_n(unsigned long iPoint) const { return Volume_n(iPoint); } /*! - * \brief Get the volume of the control volume at time n+1. + * \brief Get the volume of the control volume at time n-1. * \param[in] iPoint - Index of the point. - * \return Volume of the control volume at time n+1 + * \return Volume of the control volume at time n-1 */ inline su2double GetVolume_nM1(unsigned long iPoint) const { return Volume_nM1(iPoint); } @@ -536,10 +539,35 @@ class CPoint { void SetVolume_n(); /*! - * \brief Set the volume of the control volume at time n+1. + * \brief Set the volume of the control volume at time n-1. */ void SetVolume_nM1(); + /*! + * \brief Set the volume of the control volume at time n using n-1. + */ + void SetVolume_n_from_OldnM1(); + + /*! + * \brief Set the volume of the control volume at current time using time n. + */ + void SetVolume_from_Oldn(); + + /*! + * \brief Set the Volume to Volume_Old. + */ + void SetVolume_Old(); + + /*! + * \brief Set the Volume_n to Volume_n_Old. + */ + void SetVolume_n_Old(); + + /*! + * \brief Set the Volume_nM1 to Volume_nM1_Old. + */ + void SetVolume_nM1_Old(); + /*! * \brief Set the parent control volume of an agglomerated control volume. * \param[in] iPoint - Index of the point. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 461df39ee198..29ab797c2672 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2214,7 +2214,8 @@ enum MPI_QUANTITIES { MESH_DISPLACEMENTS = 27, /*!< \brief Mesh displacements at the interface. */ SOLUTION_TIME_N = 28, /*!< \brief Solution at time n. */ SOLUTION_TIME_N1 = 29, /*!< \brief Solution at time n-1. */ - PRIMITIVE = 30 /*!< \brief Primitive solution communication. */ + PRIMITIVE = 30, /*!< \brief Primitive solution communication. */ + SOLUTION_VEL_PRED = 31 }; /*! diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 41b60c4dd0bd..01ba04597da5 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -4834,12 +4834,6 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ Restart_Flow = false; - if (GetKind_GridMovement() != RIGID_MOTION && - GetKind_GridMovement() != NO_MOVEMENT) { - SU2_MPI::Error(string("Dynamic mesh movement currently only supported for the discrete adjoint solver for\n") + - string("GRID_MOVEMENT = RIGID_MOTION."), CURRENT_FUNCTION); - } - if (Unst_AdjointIter- long(nTimeIter) < 0){ SU2_MPI::Error(string("Invalid iteration number requested for unsteady adjoint.\n" ) + string("Make sure EXT_ITER is larger or equal than UNST_ADJOINT_ITER."), @@ -4847,7 +4841,8 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_ } /*--- If the averaging interval is not set, we average over all time-steps ---*/ - + } + if (Time_Domain) { if (Iter_Avg_Objective == 0.0) { Iter_Avg_Objective = nTimeIter; } @@ -9811,10 +9806,12 @@ void CConfig::SetMultizone(CConfig *driver_config, CConfig **config_container){ switch (config_container[iZone]->GetKind_Solver()) { case EULER: case NAVIER_STOKES: case RANS: case INC_EULER: case INC_NAVIER_STOKES: case INC_RANS: - case NEMO_EULER: case NEMO_NAVIER_STOKES: + case NEMO_EULER: case NEMO_NAVIER_STOKES: + case DISC_ADJ_EULER: case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_RANS: + case DISC_ADJ_INC_EULER: case DISC_ADJ_INC_NAVIER_STOKES: case DISC_ADJ_INC_RANS: fluid_zone = true; break; - case FEM_ELASTICITY: + case FEM_ELASTICITY: case DISC_ADJ_FEM: structural_zone = true; Relaxation = true; break; diff --git a/Common/src/geometry/dual_grid/CPoint.cpp b/Common/src/geometry/dual_grid/CPoint.cpp index d32c85c5cd00..3fda25d53222 100644 --- a/Common/src/geometry/dual_grid/CPoint.cpp +++ b/Common/src/geometry/dual_grid/CPoint.cpp @@ -68,6 +68,12 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig *config) { Volume_nM1.resize(npoint) = su2double(0.0); } + if (config->GetDynamic_Grid() && config->GetDiscrete_Adjoint()) { + Volume_Old.resize(npoint) = su2double(0.0); + Volume_n_Old.resize(npoint) = su2double(0.0); + Volume_nM1_Old.resize(npoint) = su2double(0.0); + } + if(config->GetAD_Mode() && config->GetMultizone_Problem()) { AD_InputIndex.resize(npoint,nDim) = 0; AD_OutputIndex.resize(npoint,nDim) = 0; @@ -163,6 +169,31 @@ void CPoint::SetVolume_nM1() { parallelCopy(Volume_n.size(), Volume_n.data(), Volume_nM1.data()); } +void CPoint::SetVolume_Old() { + assert(Volume_Old.size() == Volume.size()); + parallelCopy(Volume.size(), Volume.data(), Volume_Old.data()); +} + +void CPoint::SetVolume_n_Old() { + assert(Volume_n_Old.size() == Volume_n.size()); + parallelCopy(Volume_n.size(), Volume_n.data(), Volume_n_Old.data()); +} + +void CPoint::SetVolume_nM1_Old() { + assert(Volume_nM1_Old.size() == Volume_nM1.size()); + parallelCopy(Volume_nM1.size(), Volume_nM1.data(), Volume_nM1_Old.data()); +} + +void CPoint::SetVolume_n_from_OldnM1() { + assert(Volume_n.size() == Volume_nM1_Old.size()); + parallelCopy(Volume_nM1_Old.size(), Volume_nM1_Old.data(), Volume_n.data()); +} + +void CPoint::SetVolume_from_Oldn() { + assert(Volume.size() == Volume_n_Old.size()); + parallelCopy(Volume_n_Old.size(), Volume_n_Old.data(), Volume.data()); +} + void CPoint::SetCoord_n() { assert(Coord_n.size() == Coord.size()); parallelCopy(Coord.size(), Coord.data(), Coord_n.data()); diff --git a/Common/src/grid_movement/CVolumetricMovement.cpp b/Common/src/grid_movement/CVolumetricMovement.cpp index 68b35eb613a1..11f61adf35ec 100644 --- a/Common/src/grid_movement/CVolumetricMovement.cpp +++ b/Common/src/grid_movement/CVolumetricMovement.cpp @@ -438,8 +438,10 @@ void CVolumetricMovement::ComputeSolid_Wall_Distance(CGeometry *geometry, CConfi for(iMarker=0; iMarkerGetnMarker_All(); ++iMarker) { if( (config->GetMarker_All_KindBC(iMarker) == EULER_WALL || config->GetMarker_All_KindBC(iMarker) == HEAT_FLUX) || - (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) || - (config->GetMarker_All_KindBC(iMarker) == CHT_WALL_INTERFACE)) { + (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) || + (config->GetMarker_All_KindBC(iMarker) == CHT_WALL_INTERFACE) || + (config->GetMarker_All_KindBC(iMarker) == CLAMPED_BOUNDARY ) || + (config->GetMarker_All_KindBC(iMarker) == LOAD_BOUNDARY )) { nVertex_SolidWall += geometry->GetnVertex(iMarker); } } @@ -457,8 +459,10 @@ void CVolumetricMovement::ComputeSolid_Wall_Distance(CGeometry *geometry, CConfi for (iMarker=0; iMarkerGetnMarker_All(); ++iMarker) { if ( (config->GetMarker_All_KindBC(iMarker) == EULER_WALL || config->GetMarker_All_KindBC(iMarker) == HEAT_FLUX) || - (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) || - (config->GetMarker_All_KindBC(iMarker) == CHT_WALL_INTERFACE)) { + (config->GetMarker_All_KindBC(iMarker) == ISOTHERMAL) || + (config->GetMarker_All_KindBC(iMarker) == CHT_WALL_INTERFACE) || + (config->GetMarker_All_KindBC(iMarker) == CLAMPED_BOUNDARY ) || + (config->GetMarker_All_KindBC(iMarker) == LOAD_BOUNDARY )) { for (iVertex=0; iVertexGetnVertex(iMarker); ++iVertex) { iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); PointIDs[jj++] = iPoint; diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index 2d8d931776a8..25725077edc0 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -818,7 +818,8 @@ class CDriver { */ void SetSourceTerm_DispAdjoint(unsigned short iMarker, unsigned long iVertex, passivedouble val_AdjointX, passivedouble val_AdjointY, passivedouble val_AdjointZ); - + void SetSourceTerm_VelAdjoint(unsigned short iMarker, unsigned long iVertex, passivedouble val_AdjointX, + passivedouble val_AdjointY, passivedouble val_AdjointZ); /*! * \brief Get the undeformed mesh coordinates * \param[in] iMarker - Marker identifier. diff --git a/SU2_CFD/include/iteration/CAdjFluidIteration.hpp b/SU2_CFD/include/iteration/CAdjFluidIteration.hpp index 68c99d7d0dfc..0623e6437682 100644 --- a/SU2_CFD/include/iteration/CAdjFluidIteration.hpp +++ b/SU2_CFD/include/iteration/CAdjFluidIteration.hpp @@ -67,7 +67,7 @@ class CAdjFluidIteration : public CFluidIteration { void Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst) override; + unsigned short val_iInst, bool CrossTerm) override; /*! * \brief Updates the containers for the adjoint fluid system. diff --git a/SU2_CFD/include/iteration/CDiscAdjFEAIteration.hpp b/SU2_CFD/include/iteration/CDiscAdjFEAIteration.hpp index e2aac672de11..efe486c3d23f 100644 --- a/SU2_CFD/include/iteration/CDiscAdjFEAIteration.hpp +++ b/SU2_CFD/include/iteration/CDiscAdjFEAIteration.hpp @@ -90,7 +90,7 @@ class CDiscAdjFEAIteration : public CIteration { void Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst) override; + unsigned short val_iInst, bool CrossTerm) override; /*! * \brief Updates the containers for the discrete adjoint mean flow system. diff --git a/SU2_CFD/include/iteration/CDiscAdjFluidIteration.hpp b/SU2_CFD/include/iteration/CDiscAdjFluidIteration.hpp index 4545b67b461b..82f92e80081d 100644 --- a/SU2_CFD/include/iteration/CDiscAdjFluidIteration.hpp +++ b/SU2_CFD/include/iteration/CDiscAdjFluidIteration.hpp @@ -86,7 +86,7 @@ class CDiscAdjFluidIteration : public CIteration { void Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst) override; + unsigned short val_iInst, bool CrossTerm) override; /*! * \brief Updates the containers for the discrete adjoint fluid system. diff --git a/SU2_CFD/include/iteration/CDiscAdjHeatIteration.hpp b/SU2_CFD/include/iteration/CDiscAdjHeatIteration.hpp index 207c9b7f2d9e..f37c69a57afc 100644 --- a/SU2_CFD/include/iteration/CDiscAdjHeatIteration.hpp +++ b/SU2_CFD/include/iteration/CDiscAdjHeatIteration.hpp @@ -79,7 +79,7 @@ class CDiscAdjHeatIteration : public CIteration { void Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst) override; + unsigned short val_iInst, bool CrossTerm) override; /*! * \brief Perform a single iteration of the adjoint fluid system. diff --git a/SU2_CFD/include/iteration/CFEAIteration.hpp b/SU2_CFD/include/iteration/CFEAIteration.hpp index 1194d8fc3e19..1cc1e9554592 100644 --- a/SU2_CFD/include/iteration/CFEAIteration.hpp +++ b/SU2_CFD/include/iteration/CFEAIteration.hpp @@ -60,7 +60,7 @@ class CFEAIteration : public CIteration { void Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst) override; + unsigned short val_iInst, bool CrossTerm) override; /*! * \brief Iterate the structural system for a number of Inner_Iter iterations. diff --git a/SU2_CFD/include/iteration/CFEMFluidIteration.hpp b/SU2_CFD/include/iteration/CFEMFluidIteration.hpp index ae477336e56d..0266cbc5e162 100644 --- a/SU2_CFD/include/iteration/CFEMFluidIteration.hpp +++ b/SU2_CFD/include/iteration/CFEMFluidIteration.hpp @@ -75,7 +75,7 @@ class CFEMFluidIteration : public CFluidIteration { void Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst) override; + unsigned short val_iInst, bool CrossTerm) override; /*! * \brief Updates the containers for the finite element flow system. diff --git a/SU2_CFD/include/iteration/CFluidIteration.hpp b/SU2_CFD/include/iteration/CFluidIteration.hpp index dc801b1e32ca..f8d712b97366 100644 --- a/SU2_CFD/include/iteration/CFluidIteration.hpp +++ b/SU2_CFD/include/iteration/CFluidIteration.hpp @@ -67,7 +67,7 @@ class CFluidIteration : public CIteration { void Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst) override; + unsigned short val_iInst, bool CrossTerm) override; /*! * \brief Iterate the fluid system for a number of Inner_Iter iterations. diff --git a/SU2_CFD/include/iteration/CHeatIteration.hpp b/SU2_CFD/include/iteration/CHeatIteration.hpp index 6860556e526d..e3aea0be66f8 100644 --- a/SU2_CFD/include/iteration/CHeatIteration.hpp +++ b/SU2_CFD/include/iteration/CHeatIteration.hpp @@ -58,7 +58,7 @@ class CHeatIteration : public CFluidIteration { void Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst) override; + unsigned short val_iInst, bool CrossTerm) override; /*! * \brief Updates the containers for the heat system. diff --git a/SU2_CFD/include/iteration/CIteration.hpp b/SU2_CFD/include/iteration/CIteration.hpp index e20846277db9..c0107d7ca20b 100644 --- a/SU2_CFD/include/iteration/CIteration.hpp +++ b/SU2_CFD/include/iteration/CIteration.hpp @@ -136,7 +136,7 @@ class CIteration { virtual void Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst){} + unsigned short val_iInst, bool CrossTerm){} /*! * \brief A virtual member. diff --git a/SU2_CFD/include/solvers/CDiscAdjSolver.hpp b/SU2_CFD/include/solvers/CDiscAdjSolver.hpp index 52244aafd6cd..45e8f2ef2851 100644 --- a/SU2_CFD/include/solvers/CDiscAdjSolver.hpp +++ b/SU2_CFD/include/solvers/CDiscAdjSolver.hpp @@ -129,8 +129,9 @@ class CDiscAdjSolver final : public CSolver { * after tape has been evaluated. * \param[in] geometry - The geometrical definition of the problem. * \param[in] config - The particular config. + * \param[in] CrossTerm - Boolean to determine if this is a cross term extraction. */ - void ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config) override; + void ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm) override; /*! * \brief A virtual member. diff --git a/SU2_CFD/include/solvers/CFEASolver.hpp b/SU2_CFD/include/solvers/CFEASolver.hpp index 6602d21ccbba..993e284ddeac 100644 --- a/SU2_CFD/include/solvers/CFEASolver.hpp +++ b/SU2_CFD/include/solvers/CFEASolver.hpp @@ -449,7 +449,10 @@ class CFEASolver : public CSolver { CNumerics *numerics, const CConfig *config, unsigned short val_marker) final; - + void BC_Velocity(CGeometry *geometry, + CNumerics *numerics, + const CConfig *config, + unsigned short val_marker) final; /*! * \brief Iterate using an implicit Newmark solver. * \param[in] geometry - Geometrical definition of the problem. @@ -638,6 +641,7 @@ class CFEASolver : public CSolver { * \param[in] config - Configuration of the problem. */ void PredictStruct_Displacement(CGeometry *geometry, CConfig *config) final; + void PredictStruct_Velocity(CGeometry *geometry, CConfig *config) final; /*! * \brief Computation of Aitken's coefficient. diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 26b4bb7850b7..9d542f011b09 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -1369,7 +1369,7 @@ void CFVMFlowSolverBase::Pressure_Forces(const CGeometry* geometr RefTemp = Temperature_Inf; RefDensity = Density_Inf; - if (dynamic_grid) { + if (dynamic_grid && !config->GetMultizone_Problem()) { Mach2Vel = sqrt(Gamma * Gas_Constant * RefTemp); Mach_Motion = config->GetMach_Motion(); RefVel2 = (Mach_Motion * Mach2Vel) * (Mach_Motion * Mach2Vel); diff --git a/SU2_CFD/include/solvers/CMeshSolver.hpp b/SU2_CFD/include/solvers/CMeshSolver.hpp index e94b4e1a222a..95b2e2ccd0af 100644 --- a/SU2_CFD/include/solvers/CMeshSolver.hpp +++ b/SU2_CFD/include/solvers/CMeshSolver.hpp @@ -87,6 +87,7 @@ class CMeshSolver final : public CFEASolver { * \param[in] config - Definition of the particular problem. */ void ComputeGridVelocity(CGeometry *geometry, CConfig *config); + void ComputeGridVelocity_FromBoundary(CGeometry **geometry, CNumerics **numerics, CConfig *config); /*! * \brief Update the coarse multigrid levels after the grid movement. @@ -101,6 +102,7 @@ class CMeshSolver final : public CFEASolver { * \param[in] config - Definition of the particular problem. */ void SetBoundaryDisplacements(CGeometry *geometry, CNumerics *numerics, CConfig *config); + void SetBoundaryVelocities(CGeometry *geometry, CNumerics *numerics, CConfig *config); public: /*! diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index fd37e55262b5..d80be5386148 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -985,6 +985,10 @@ class CSolver { CNumerics *numerics, const CConfig *config, unsigned short val_marker) { } + inline virtual void BC_Velocity(CGeometry *geometry, + CNumerics *numerics, + const CConfig *config, + unsigned short val_marker) { } /*! * \brief A virtual member. @@ -3514,6 +3518,8 @@ class CSolver { */ inline virtual void PredictStruct_Displacement(CGeometry *geometry, CConfig *config) { } + inline virtual void PredictStruct_Velocity(CGeometry *geometry, + CConfig *config) { } /*! * \brief A virtual member. @@ -3775,6 +3781,15 @@ class CSolver { * \param[in] config - The particular config. */ inline virtual void ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config){} + + /*! + * \brief A virtual member. + * \param[in] geometry - The geometrical definition of the problem. + * \param[in] solver_container - The solver container holding all solutions. + * \param[in] config - The particular config. + * \param[in] CrossTerm - Boolean to determine if this is a cross term extraction. + */ + inline virtual void ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm){} /*! * \brief A virtual member. diff --git a/SU2_CFD/include/variables/CDiscAdjFEABoundVariable.hpp b/SU2_CFD/include/variables/CDiscAdjFEABoundVariable.hpp index ede224e14b6a..4938c9372349 100644 --- a/SU2_CFD/include/variables/CDiscAdjFEABoundVariable.hpp +++ b/SU2_CFD/include/variables/CDiscAdjFEABoundVariable.hpp @@ -43,6 +43,7 @@ class CDiscAdjFEABoundVariable final : public CDiscAdjFEAVariable { MatrixType FlowTraction_Sens; /*!< \brief Adjoint of the flow tractions. */ MatrixType SourceTerm_DispAdjoint; /*!< \brief Source term applied into the displacement adjoint coming from external solvers. */ + MatrixType SourceTerm_VelAdjoint; CVertexMap VertexMap; /*!< \brief Object that controls accesses to the variables of this class. */ @@ -101,6 +102,10 @@ class CDiscAdjFEABoundVariable final : public CDiscAdjFEAVariable { if (!VertexMap.GetVertexIndex(iPoint)) return; SourceTerm_DispAdjoint(iPoint,iDim) = val; } + inline void SetSourceTerm_VelAdjoint(unsigned long iPoint, unsigned long iDim, su2double val) override { + if (!VertexMap.GetVertexIndex(iPoint)) return; + SourceTerm_VelAdjoint(iPoint,iDim) = val; + } /*! * \brief Get the source term applied into the displacement adjoint coming from external solvers @@ -111,6 +116,10 @@ class CDiscAdjFEABoundVariable final : public CDiscAdjFEAVariable { if (!VertexMap.GetVertexIndex(iPoint)) return 0.0; return SourceTerm_DispAdjoint(iPoint,iDim); } + inline su2double GetSourceTerm_VelAdjoint(unsigned long iPoint, unsigned long iDim) const override { + if (!VertexMap.GetVertexIndex(iPoint)) return 0.0; + return SourceTerm_VelAdjoint(iPoint,iDim); + } /*! * \brief Get whether a node is on the boundary diff --git a/SU2_CFD/include/variables/CDiscAdjFEAVariable.hpp b/SU2_CFD/include/variables/CDiscAdjFEAVariable.hpp index 289bd1dc7873..209e741d820d 100644 --- a/SU2_CFD/include/variables/CDiscAdjFEAVariable.hpp +++ b/SU2_CFD/include/variables/CDiscAdjFEAVariable.hpp @@ -39,6 +39,7 @@ class CDiscAdjFEAVariable : public CVariable { protected: MatrixType Sensitivity; /* Vector holding the derivative of target functional with respect to the coordinates at this node*/ + MatrixType Sensitivity_Old; /* Previous time sensitivity holder since inner iterations in FSI problems overwrite sensitivity*/ MatrixType Solution_Direct; MatrixType Dynamic_Derivative; @@ -87,6 +88,8 @@ class CDiscAdjFEAVariable : public CVariable { */ inline void SetSensitivity(unsigned long iPoint, unsigned long iDim, su2double val) final { Sensitivity(iPoint,iDim) = val; } + inline void SetSensitivity_Old(unsigned long iPoint, unsigned long iDim, su2double val) final { Sensitivity_Old(iPoint,iDim) = val; } + /*! * \brief Get the Sensitivity at the node * \param[in] iDim - spacial component @@ -94,6 +97,8 @@ class CDiscAdjFEAVariable : public CVariable { */ inline su2double GetSensitivity(unsigned long iPoint, unsigned long iDim) const final { return Sensitivity(iPoint,iDim);} + inline su2double GetSensitivity_Old(unsigned long iPoint, unsigned long iDim) const final { return Sensitivity_Old(iPoint,iDim);} + inline void SetDynamic_Derivative(unsigned long iPoint, unsigned long iVar, su2double der) final { Dynamic_Derivative(iPoint,iVar) = der; } diff --git a/SU2_CFD/include/variables/CDiscAdjVariable.hpp b/SU2_CFD/include/variables/CDiscAdjVariable.hpp index aa7185b6f753..562f86422a14 100644 --- a/SU2_CFD/include/variables/CDiscAdjVariable.hpp +++ b/SU2_CFD/include/variables/CDiscAdjVariable.hpp @@ -38,6 +38,7 @@ class CDiscAdjVariable final : public CVariable { private: MatrixType Sensitivity; /* Vector holding the derivative of target functional with respect to the coordinates at this node*/ + MatrixType Sensitivity_Old; /* Previous time sensitivity holder since inner iterations in FSI problems overwrite sensitivity*/ MatrixType Solution_Direct; MatrixType DualTime_Derivative; MatrixType DualTime_Derivative_n; @@ -71,6 +72,8 @@ class CDiscAdjVariable final : public CVariable { */ inline void SetSensitivity(unsigned long iPoint, unsigned long iDim, su2double val) override { Sensitivity(iPoint,iDim) = val;} + inline void SetSensitivity_Old(unsigned long iPoint, unsigned long iDim, su2double val) override { Sensitivity_Old(iPoint,iDim) = val;} + /*! * \brief Get the Sensitivity at the node * \param[in] iDim - spacial component @@ -78,6 +81,8 @@ class CDiscAdjVariable final : public CVariable { */ inline su2double GetSensitivity(unsigned long iPoint, unsigned long iDim) const override { return Sensitivity(iPoint,iDim); } + inline su2double GetSensitivity_Old(unsigned long iPoint, unsigned long iDim) const override { return Sensitivity_Old(iPoint,iDim); } + inline void SetDual_Time_Derivative(unsigned long iPoint, unsigned long iVar, su2double der) override { DualTime_Derivative(iPoint,iVar) = der; } inline void SetDual_Time_Derivative_n(unsigned long iPoint, unsigned long iVar, su2double der) override { DualTime_Derivative_n(iPoint,iVar) = der; } diff --git a/SU2_CFD/include/variables/CFEAVariable.hpp b/SU2_CFD/include/variables/CFEAVariable.hpp index 9c3838df25ce..0da7db745661 100644 --- a/SU2_CFD/include/variables/CFEAVariable.hpp +++ b/SU2_CFD/include/variables/CFEAVariable.hpp @@ -52,12 +52,25 @@ class CFEAVariable : public CVariable { MatrixType Solution_Accel_time_n; /*!< \brief Acceleration of the nodes at time n. */ MatrixType Solution_Pred; /*!< \brief Predictor of the solution for FSI purposes */ + MatrixType Solution_Vel_Pred; MatrixType Solution_Pred_Old; /*!< \brief Predictor of the solution at time n for FSI purposes */ MatrixType Reference_Geometry; /*!< \brief Reference solution for optimization problems */ MatrixType Prestretch; /*!< \brief Prestretch geometry */ + su2matrix AD_Vel_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ + su2matrix AD_Vel_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + + su2matrix AD_Vel_Time_n_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ + su2matrix AD_Vel_Time_n_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + + su2matrix AD_Accel_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ + su2matrix AD_Accel_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + + su2matrix AD_Accel_Time_n_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ + su2matrix AD_Accel_Time_n_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + /*! * \brief Constructor of the class. * \note This class is not supposed to be instantiated, it is only a building block for CFEABoundVariable @@ -257,6 +270,9 @@ class CFEAVariable : public CVariable { inline void SetSolution_Pred(unsigned long iPoint) final { for (unsigned long iVar = 0; iVar < nVar; iVar++) Solution_Pred(iPoint,iVar) = Solution(iPoint,iVar); } + inline void SetSolution_Vel_Pred(unsigned long iPoint) final { + for (unsigned long iVar = 0; iVar < nVar; iVar++) Solution_Vel_Pred(iPoint,iVar) = Solution_Vel(iPoint,iVar); + } /*! * \brief Set the value of the old solution. @@ -265,6 +281,9 @@ class CFEAVariable : public CVariable { inline void SetSolution_Pred(unsigned long iPoint, const su2double *val_solution_pred) final { for (unsigned long iVar = 0; iVar < nVar; iVar++) Solution_Pred(iPoint,iVar) = val_solution_pred[iVar]; } + inline void SetSolution_Vel_Pred(unsigned long iPoint, const su2double *val_solution_pred) final { + for (unsigned long iVar = 0; iVar < nVar; iVar++) Solution_Vel_Pred(iPoint,iVar) = val_solution_pred[iVar]; + } /*! * \brief Set the value of the predicted solution. @@ -274,6 +293,9 @@ class CFEAVariable : public CVariable { inline void SetSolution_Pred(unsigned long iPoint, unsigned long iVar, su2double val_solution_pred) final { Solution_Pred(iPoint,iVar) = val_solution_pred; } + inline void SetSolution_Vel_Pred(unsigned long iPoint, unsigned long iVar, su2double val_solution_pred) final { + Solution_Vel_Pred(iPoint,iVar) = val_solution_pred; + } /*! * \brief Get the value of the solution predictor. @@ -281,12 +303,14 @@ class CFEAVariable : public CVariable { * \return Pointer to the old solution vector. */ inline su2double GetSolution_Pred(unsigned long iPoint, unsigned long iVar) const final { return Solution_Pred(iPoint,iVar); } + inline su2double GetSolution_Vel_Pred(unsigned long iPoint, unsigned long iVar) const final { return Solution_Vel_Pred(iPoint,iVar); } /*! * \brief Get the solution at time n. * \return Pointer to the solution (at time n) vector. */ inline su2double *GetSolution_Pred(unsigned long iPoint) final { return Solution_Pred[iPoint]; } + inline su2double *GetSolution_Vel_Pred(unsigned long iPoint) final { return Solution_Vel_Pred[iPoint]; } /*! * \brief Set the value of the solution predictor. @@ -380,29 +404,29 @@ class CFEAVariable : public CVariable { * \brief Register the variables in the solution time_n array as input/output variable. * \param[in] input - input or output variables. */ - void Register_femSolution_time_n() final; + void Register_femSolution_time_n(bool input, bool push_index) final; /*! * \brief Register the variables in the velocity array as input/output variable. * \param[in] input - input or output variables. */ - void RegisterSolution_Vel(bool input) final; + void RegisterSolution_Vel(bool input, bool push_index) final; /*! * \brief Register the variables in the velocity time_n array as input/output variable. */ - void RegisterSolution_Vel_time_n() final; + void RegisterSolution_Vel_time_n(bool input, bool push_index) final; /*! * \brief Register the variables in the acceleration array as input/output variable. * \param[in] input - input or output variables. */ - void RegisterSolution_Accel(bool input) final; + void RegisterSolution_Accel(bool input, bool push_index) final; /*! * \brief Register the variables in the acceleration time_n array as input/output variable. */ - void RegisterSolution_Accel_time_n() final; + void RegisterSolution_Accel_time_n(bool input, bool push_index) final; /*! * \brief Set the velocity adjoint values of the solution. @@ -422,6 +446,11 @@ class CFEAVariable : public CVariable { adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_Vel(iPoint,iVar)); } + inline void GetAdjointSolution_Vel_LocalIndex(unsigned long iPoint, su2double *adj_sol) const final { + for (unsigned long iVar = 0; iVar < nVar; iVar++) + adj_sol[iVar] = AD::GetDerivative(AD_Vel_InputIndex(iPoint,iVar)); + } + /*! * \brief Set the velocity adjoint values of the solution at time n. * \param[in] adj_sol - The adjoint values of the solution. @@ -440,6 +469,11 @@ class CFEAVariable : public CVariable { adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_Vel_time_n(iPoint,iVar)); } + inline void GetAdjointSolution_Vel_time_n_LocalIndex(unsigned long iPoint, su2double *adj_sol) const final { + for (unsigned long iVar = 0; iVar < nVar; iVar++) + adj_sol[iVar] = AD::GetDerivative(AD_Vel_Time_n_InputIndex(iPoint,iVar)); + } + /*! * \brief Set the acceleration adjoint values of the solution. * \param[in] adj_sol - The adjoint values of the solution. @@ -458,6 +492,11 @@ class CFEAVariable : public CVariable { adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_Accel(iPoint,iVar)); } + inline void GetAdjointSolution_Accel_LocalIndex(unsigned long iPoint, su2double *adj_sol) const final { + for (unsigned long iVar = 0; iVar < nVar; iVar++) + adj_sol[iVar] = AD::GetDerivative(AD_Accel_InputIndex(iPoint,iVar)); + } + /*! * \brief Set the acceleration adjoint values of the solution at time n. * \param[in] adj_sol - The adjoint values of the solution. @@ -476,4 +515,8 @@ class CFEAVariable : public CVariable { adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_Accel_time_n(iPoint,iVar)); } + inline void GetAdjointSolution_Accel_time_n_LocalIndex(unsigned long iPoint, su2double *adj_sol) const final { + for (unsigned long iVar = 0; iVar < nVar; iVar++) + adj_sol[iVar] = AD::GetDerivative(AD_Accel_Time_n_InputIndex(iPoint,iVar)); + } }; diff --git a/SU2_CFD/include/variables/CMeshBoundVariable.hpp b/SU2_CFD/include/variables/CMeshBoundVariable.hpp index aadcae288604..766f8039c466 100644 --- a/SU2_CFD/include/variables/CMeshBoundVariable.hpp +++ b/SU2_CFD/include/variables/CMeshBoundVariable.hpp @@ -35,7 +35,7 @@ class CMeshBoundVariable final : public CMeshVariable { private: MatrixType Boundary_Displacement; /*!< \brief Store the reference coordinates of the mesh. */ - + MatrixType Boundary_Velocity; CVertexMap VertexMap; /*!< \brief Object that controls accesses to the variables of this class. */ public: @@ -87,6 +87,34 @@ class CMeshBoundVariable final : public CMeshVariable { Boundary_Displacement(iPoint,iDim) = val_BoundDisp; } +/*! + * \brief Get the value of the displacement imposed at the boundary. + * \return Value of the boundary velocity. + */ + inline su2double GetBound_Vel(unsigned long iPoint, unsigned long iDim) const override { + if (!VertexMap.GetVertexIndex(iPoint)) return 0.0; + return Boundary_Velocity(iPoint,iDim); + } + + /*! + * \brief Set the boundary displacements. + * \param[in] val_BoundVel - Pointer to the boundary velocities. + */ + inline void SetBound_Vel(unsigned long iPoint, const su2double *val_BoundVel) override { + if (!VertexMap.GetVertexIndex(iPoint)) return; + for (unsigned long iDim = 0; iDim < nDim; iDim++) Boundary_Velocity(iPoint,iDim) = val_BoundVel[iDim]; + } + + /*! + * \brief Set the boundary velocity. + * \param[in] iDim - Index of the dimension of interest. + * \param[in] val_BoundVel - Value of the boundary velocities. + */ + inline void SetBound_Vel(unsigned long iPoint, unsigned long iDim, su2double val_BoundVel) override { + if (!VertexMap.GetVertexIndex(iPoint)) return; + Boundary_Velocity(iPoint,iDim) = val_BoundVel; + } + /*! * \brief Register the boundary displacements of the mesh. * \param[in] input - Defines whether we are registering the variable as input or as output. diff --git a/SU2_CFD/include/variables/CVariable.hpp b/SU2_CFD/include/variables/CVariable.hpp index b103e17f4709..6b7ee864f2c7 100644 --- a/SU2_CFD/include/variables/CVariable.hpp +++ b/SU2_CFD/include/variables/CVariable.hpp @@ -98,6 +98,12 @@ class CVariable { su2matrix AD_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ su2matrix AD_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + su2matrix AD_Time_n_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ + su2matrix AD_Time_n_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + + su2matrix AD_Time_n1_InputIndex; /*!< \brief Indices of Solution variables in the adjoint vector. */ + su2matrix AD_Time_n1_OutputIndex; /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */ + unsigned long nPoint = 0; /*!< \brief Number of points in the domain. */ unsigned long nDim = 0; /*!< \brief Number of dimension of the problem. */ unsigned long nVar = 0; /*!< \brief Number of variables of the problem. */ @@ -2244,18 +2250,21 @@ class CVariable { * \brief A virtual member. Set the value of the solution predictor. */ inline virtual void SetSolution_Pred(unsigned long iPoint) {} + inline virtual void SetSolution_Vel_Pred(unsigned long iPoint) {} /*! * \brief A virtual member. Set the value of the old solution. * \param[in] solution_pred - Pointer to the residual vector. */ inline virtual void SetSolution_Pred(unsigned long iPoint, const su2double *solution_pred) {} + inline virtual void SetSolution_Vel_Pred(unsigned long iPoint, const su2double *solution_pred) {} /*! * \brief A virtual member. Set the value of the solution predicted. * \param[in] solution_old - Pointer to the residual vector. */ inline virtual void SetSolution_Pred(unsigned long iPoint, unsigned long iVar, su2double solution_pred) {} + inline virtual void SetSolution_Vel_Pred(unsigned long iPoint, unsigned long iVar, su2double solution_pred) {} /*! * \brief A virtual member. Get the value of the solution predictor. @@ -2263,12 +2272,14 @@ class CVariable { * \return Pointer to the old solution vector. */ inline virtual su2double GetSolution_Pred(unsigned long iPoint, unsigned long iVar) const { return 0.0; } + inline virtual su2double GetSolution_Vel_Pred(unsigned long iPoint, unsigned long iVar) const { return 0.0; } /*! * \brief A virtual member. Get the solution at time n. * \return Pointer to the solution (at time n) vector. */ inline virtual su2double *GetSolution_Pred(unsigned long iPoint) {return nullptr; } + inline virtual su2double *GetSolution_Vel_Pred(unsigned long iPoint) {return nullptr; } /*! * \brief A virtual member. Set the value of the solution predictor. @@ -2379,13 +2390,14 @@ class CVariable { * \return Value of the boundary displacement. */ inline virtual su2double GetBound_Disp(unsigned long iPoint, unsigned long iDim) const { return 0.0; } + inline virtual su2double GetBound_Vel(unsigned long iPoint, unsigned long iDim) const { return 0.0; } /*! * \brief A virtual member. Set the boundary displacement. * \param[in] val_BoundDisp - Pointer to the boundary displacements. */ inline virtual void SetBound_Disp(unsigned long iPoint, const su2double *val_BoundDisp) { } - + inline virtual void SetBound_Vel(unsigned long iPoint, const su2double *val_BoundDisp) { } /*! * \brief A virtual member. Set the boundary displacement. @@ -2393,6 +2405,7 @@ class CVariable { * \param[in] val_BoundDisp - Value of the boundary displacements. */ inline virtual void SetBound_Disp(unsigned long iPoint, unsigned long iDim, const su2double val_BoundDisp) { } + inline virtual void SetBound_Vel(unsigned long iPoint, unsigned long iDim, const su2double val_BoundDisp) { } /*! * \brief A virtual member. Get the value of the displacement imposed at the boundary. @@ -2433,27 +2446,27 @@ class CVariable { /*! * \brief A virtual member. */ - inline virtual void Register_femSolution_time_n() {} + inline virtual void Register_femSolution_time_n(bool input, bool push_index) {} /*! * \brief A virtual member. */ - inline virtual void RegisterSolution_Vel(bool input) {} + inline virtual void RegisterSolution_Vel(bool input, bool push_index) {} /*! * \brief A virtual member. */ - inline virtual void RegisterSolution_Vel_time_n() {} + inline virtual void RegisterSolution_Vel_time_n(bool input, bool push_index) {} /*! * \brief A virtual member. */ - inline virtual void RegisterSolution_Accel(bool input) {} + inline virtual void RegisterSolution_Accel(bool input, bool push_index) {} /*! * \brief A virtual member. */ - inline virtual void RegisterSolution_Accel_time_n() {} + inline virtual void RegisterSolution_Accel_time_n(bool input, bool push_index) {} /*! * \brief A virtual member. @@ -2474,7 +2487,7 @@ class CVariable { * \brief A virtual member. */ inline virtual void GetAdjointSolution_Vel(unsigned long iPoint, su2double *adj_sol) const {} - + inline virtual void GetAdjointSolution_Vel_LocalIndex(unsigned long iPoint, su2double *adj_sol) const {} /*! * \brief A virtual member. */ @@ -2484,6 +2497,7 @@ class CVariable { * \brief A virtual member. */ inline virtual void GetAdjointSolution_Vel_time_n(unsigned long iPoint, su2double *adj_sol) const {} + inline virtual void GetAdjointSolution_Vel_time_n_LocalIndex(unsigned long iPoint, su2double *adj_sol) const {} /*! * \brief A virtual member. @@ -2494,6 +2508,7 @@ class CVariable { * \brief A virtual member. */ inline virtual void GetAdjointSolution_Accel(unsigned long iPoint, su2double *adj_sol) const {} + inline virtual void GetAdjointSolution_Accel_LocalIndex(unsigned long iPoint, su2double *adj_sol) const {} /*! * \brief A virtual member. @@ -2504,6 +2519,7 @@ class CVariable { * \brief A virtual member. */ inline virtual void GetAdjointSolution_Accel_time_n(unsigned long iPoint, su2double *adj_sol) const {} + inline virtual void GetAdjointSolution_Accel_time_n_LocalIndex(unsigned long iPoint, su2double *adj_sol) const {} /*! * \brief Register the variables in the solution array as input/output variable. @@ -2515,12 +2531,12 @@ class CVariable { /*! * \brief Register the variables in the solution_time_n array as input/output variable. */ - void RegisterSolution_time_n(); + void RegisterSolution_time_n(bool push_index); /*! * \brief Register the variables in the solution_time_n1 array as input/output variable. */ - void RegisterSolution_time_n1(); + void RegisterSolution_time_n1(bool push_index); /*! * \brief Set the adjoint values of the solution. @@ -2576,6 +2592,11 @@ class CVariable { adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_time_n(iPoint,iVar)); } + inline void GetAdjointSolution_time_n_LocalIndex(unsigned long iPoint, su2double *adj_sol) const { + for (unsigned long iVar = 0; iVar < nVar; iVar++) + adj_sol[iVar] = AD::GetDerivative(AD_Time_n_InputIndex(iPoint,iVar)); + } + /*! * \brief Set the adjoint values of the solution at time n-1. * \param[in] adj_sol - The adjoint values of the solution. @@ -2594,6 +2615,11 @@ class CVariable { adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_time_n1(iPoint,iVar)); } + inline void GetAdjointSolution_time_n1_LocalIndex(unsigned long iPoint, su2double *adj_sol) const { + for (unsigned long iVar = 0; iVar < nVar; iVar++) + adj_sol[iVar] = AD::GetDerivative(AD_Time_n1_InputIndex(iPoint,iVar)); + } + /*! * \brief Set the sensitivity at the node * \param[in] iDim - spacial component @@ -2601,6 +2627,8 @@ class CVariable { */ inline virtual void SetSensitivity(unsigned long iPoint, unsigned long iDim, su2double val) {} + inline virtual void SetSensitivity_Old(unsigned long iPoint, unsigned long iDim, su2double val) {} + /*! * \brief Get the Sensitivity at the node * \param[in] iDim - spacial component @@ -2608,6 +2636,8 @@ class CVariable { */ inline virtual su2double GetSensitivity(unsigned long iPoint, unsigned long iDim) const { return 0.0; } + inline virtual su2double GetSensitivity_Old(unsigned long iPoint, unsigned long iDim) const { return 0.0; } + inline virtual void SetDual_Time_Derivative(unsigned long iPoint, unsigned long iVar, su2double der) {} inline virtual void SetDual_Time_Derivative_n(unsigned long iPoint, unsigned long iVar, su2double der) {} @@ -2691,6 +2721,7 @@ class CVariable { * \param[in] val - value of the source term */ virtual void SetSourceTerm_DispAdjoint(unsigned long iPoint, unsigned long iDim, su2double val) { } + virtual void SetSourceTerm_VelAdjoint(unsigned long iPoint, unsigned long iDim, su2double val) { } /*! * \brief Get the source term applied into the displacement adjoint coming from external solvers @@ -2698,5 +2729,6 @@ class CVariable { * \return value of the source term */ virtual su2double GetSourceTerm_DispAdjoint(unsigned long iPoint, unsigned long iDim) const { return 0.0; } + virtual su2double GetSourceTerm_VelAdjoint(unsigned long iPoint, unsigned long iDim) const { return 0.0; } }; diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 61516cce31c2..ae2120106a33 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -122,6 +122,8 @@ CDiscAdjMultizoneDriver::~CDiscAdjMultizoneDriver(){ void CDiscAdjMultizoneDriver::StartSolver() { + bool time_domain = driver_config->GetTime_Domain(); + /*--- Main external loop of the solver. Runs for the number of time steps required. ---*/ if (rank == MASTER_NODE) { @@ -129,42 +131,61 @@ void CDiscAdjMultizoneDriver::StartSolver() { cout << "\nSimulation Run using the Discrete Adjoint Multizone Driver" << endl; - if (driver_config->GetTime_Domain()) - SU2_MPI::Error("The discrete adjoint multizone driver is not ready for unsteady computations yet.", - CURRENT_FUNCTION); + if (time_domain) + cout << "The simulation will run for " << driver_config->GetnTime_Iter() << " time steps." << endl; } - for (iZone = 0; iZone < nZone; iZone++){ + while ( TimeIter < driver_config->GetnTime_Iter()) { + + for (iZone = 0; iZone < nZone; iZone++) { - /*--- Set the value of the external iteration to TimeIter. -------------------------------------*/ - /*--- TODO: This should be generalised for an homogeneous criteria throughout the code. --------*/ - config_container[iZone]->SetTimeIter(0); + /*--- Set current time iteration ---*/ + config_container[iZone]->SetTimeIter(TimeIter); - } + if (time_domain) + config_container[iZone]->SetPhysicalTime(static_cast(TimeIter)*config_container[iZone]->GetDelta_UnstTimeND()); + else + config_container[iZone]->SetPhysicalTime(0.0); + } - /*--- Size and initialize the matrix of cross-terms. ---*/ + /*--- Size and initialize the matrix of cross-terms. ---*/ - InitializeCrossTerms(); + InitializeCrossTerms(); - /*--- We directly start the (steady-state) discrete adjoint computation. ---*/ + /*--- We directly start the discrete adjoint computation. ---*/ - Run(); + Run(); - /*--- Output the solution in files. ---*/ + /*--- Output the solution in files for each time iteration. ---*/ - Output(TimeIter); + Output(TimeIter); + TimeIter++; + } } void CDiscAdjMultizoneDriver::Run() { unsigned long wrt_sol_freq = 9999; unsigned long nOuterIter = driver_config->GetnOuter_Iter(); + bool time_domain = driver_config->GetTime_Domain(); vector > fixPtCorrector(nZone); + // Reset external and solution for (iZone = 0; iZone < nZone; iZone++) { + for (unsigned short iSol=0; iSol < MAX_SOLS; iSol++) { + auto solver = solver_container[iZone][INST_0][MESH_0][iSol]; + if (solver != nullptr) { + solver->GetNodes()->SetExternalZero(); + } + } + Set_Solution_To_BGSSolution_k(iZone); + } - wrt_sol_freq = min(wrt_sol_freq, config_container[iZone]->GetVolume_Wrt_Freq()); + for (iZone = 0; iZone < nZone; iZone++) { + + if (!time_domain) + wrt_sol_freq = min(wrt_sol_freq, config_container[iZone]->GetVolume_Wrt_Freq()); iteration_container[iZone][INST_0]->Preprocess(output_container[iZone], integration_container, geometry_container, solver_container, numerics_container, config_container, surface_movement, @@ -203,7 +224,7 @@ void CDiscAdjMultizoneDriver::Run() { iteration_container[iZone][INST_0]->Iterate(output_container[iZone], integration_container, geometry_container, solver_container, numerics_container, config_container, - surface_movement, grid_movement, FFDBox, iZone, INST_0); + surface_movement, grid_movement, FFDBox, iZone, INST_0, false); Add_Solution_To_External(iZone); for (unsigned short iSol=0; iSol < MAX_SOLS; iSol++) { @@ -303,10 +324,9 @@ void CDiscAdjMultizoneDriver::Run() { ComputeAdjoints(iZone, eval_transfer); /*--- Extracting adjoints for solvers in iZone w.r.t. to outputs in iZone (diagonal part). ---*/ - iteration_container[iZone][INST_0]->Iterate(output_container[iZone], integration_container, geometry_container, solver_container, numerics_container, config_container, - surface_movement, grid_movement, FFDBox, iZone, INST_0); + surface_movement, grid_movement, FFDBox, iZone, INST_0, false); /*--- Use QN driver to improve the solution. ---*/ @@ -343,7 +363,7 @@ void CDiscAdjMultizoneDriver::Run() { config_container[jZone]->SetInnerIter(0); iteration_container[jZone][INST_0]->Iterate(output_container[jZone], integration_container, geometry_container, solver_container, numerics_container, config_container, - surface_movement, grid_movement, FFDBox, jZone, INST_0); + surface_movement, grid_movement, FFDBox, jZone, INST_0, true); /*--- Extract the cross-term performing a relaxed update of it and of the sum (External) for jZone. ---*/ @@ -359,7 +379,7 @@ void CDiscAdjMultizoneDriver::Run() { /*--- Set the multizone output. ---*/ - driver_output->SetMultizoneHistory_Output(output_container, config_container, driver_config, 0, iOuterIter); + driver_output->SetMultizoneHistory_Output(output_container, config_container, driver_config, TimeIter, iOuterIter); /*--- Check for convergence. ---*/ @@ -376,6 +396,46 @@ void CDiscAdjMultizoneDriver::Run() { if (checkSensitivity) EvaluateSensitivities(iOuterIter, StopCalc); } + + /*--- Add current time sensitivity ---*/ + /*--- Cumulative sensitivity stored here in order to exclude summation within outer iteration ---*/ + for (iZone = 0; iZone < nZone; iZone++) { + + auto solvers = solver_container[iZone][INST_0][MESH_0]; + auto geometry = geometry_container[iZone][INST_0][MESH_0]; + int IDX_SOL;//unsigned short or int? + + switch (config_container[iZone]->GetKind_Solver()) { + case DISC_ADJ_EULER: case DISC_ADJ_NAVIER_STOKES: case DISC_ADJ_RANS: + case DISC_ADJ_INC_EULER: case DISC_ADJ_INC_NAVIER_STOKES: case DISC_ADJ_INC_RANS: + if(config_container[iZone]->GetDeform_Mesh()) IDX_SOL = ADJMESH_SOL; + else IDX_SOL = ADJFLOW_SOL; + break; + case DISC_ADJ_HEAT: + IDX_SOL = ADJHEAT_SOL; + break; + case DISC_ADJ_FEM: + IDX_SOL = ADJFEA_SOL; + break; + } + + su2double Sensitivity; + for (unsigned long iPoint = 0; iPoint < geometry->GetnPoint(); iPoint++) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + + /*--- Current time iteration sensitivity ---*/ + Sensitivity = solvers[IDX_SOL]->GetNodes()->GetSensitivity(iPoint, iDim); + + /*--- Update old sensitivity container ---*/ + solvers[IDX_SOL]->GetNodes()->SetSensitivity_Old(iPoint, iDim, Sensitivity + + solvers[IDX_SOL]->GetNodes()->GetSensitivity_Old(iPoint, iDim)); + + /*--- Update sensitivity ---*/ + solvers[IDX_SOL]->GetNodes()->SetSensitivity( + iPoint, iDim, solvers[IDX_SOL]->GetNodes()->GetSensitivity_Old(iPoint, iDim)); + } + } + } } void CDiscAdjMultizoneDriver::EvaluateSensitivities(unsigned long iOuterIter, bool StopCalc) { @@ -455,11 +515,13 @@ void CDiscAdjMultizoneDriver::EvaluateSensitivities(unsigned long iOuterIter, bo AD::ClearAdjoints(); - for (iZone = 0; iZone < nZone; iZone++) { + if (!driver_config->GetTime_Domain()) {//Output files for each time iteration are handled in Output(TimeIter) + for (iZone = 0; iZone < nZone; iZone++) { - output_container[iZone]->SetResult_Files(geometry_container[iZone][INST_0][MESH_0], - config_container[iZone], - solver_container[iZone][INST_0][MESH_0], iOuterIter, StopCalc); + output_container[iZone]->SetResult_Files(geometry_container[iZone][INST_0][MESH_0], + config_container[iZone], + solver_container[iZone][INST_0][MESH_0], iOuterIter, StopCalc); + } } } @@ -505,7 +567,7 @@ void CDiscAdjMultizoneDriver::SetRecording(unsigned short kind_recording, Kind_T unsigned short type_recording = kind_recording; - if (Has_Deformation(iZone) && (kind_recording == MESH_COORDS)) { + if (Has_Deformation(iZone) && (kind_recording == MESH_COORDS) && config_container[iZone]->GetDeform_Mesh()) { type_recording = MESH_DEFORM; } @@ -595,7 +657,7 @@ void CDiscAdjMultizoneDriver::DirectIteration(unsigned short iZone, unsigned sho /*--- Iterate the zone as a block a single time ---*/ direct_iteration[iZone][INST_0]->Iterate(output_container[iZone], integration_container, geometry_container, solver_container, numerics_container, config_container, - surface_movement, grid_movement, FFDBox, iZone, INST_0); + surface_movement, grid_movement, FFDBox, iZone, INST_0, false); /*--- Print residuals in the first iteration ---*/ @@ -932,7 +994,7 @@ void CDiscAdjMultizoneDriver::HandleDataTransfer() { for(iZone = 0; iZone < nZone; iZone++) { /*--- In principle, the mesh does not need to be updated ---*/ - bool DeformMesh = false; + bool DeformMesh = (false || (config_container[iZone]->GetDeform_Mesh() && config_container[iZone]->GetSurface_Movement(DEFORMING))); /*--- Transfer from all the remaining zones ---*/ for (unsigned short jZone = 0; jZone < nZone; jZone++){ diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp index b1561b048f94..0da9872a7cc9 100644 --- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp @@ -144,7 +144,7 @@ void CDiscAdjSinglezoneDriver::Preprocess(unsigned long TimeIter) { *--- respect to the conservative variables. Since these derivatives do not change in the steady state case *--- we only have to record if the current recording is different from the main variables. ---*/ - if (RecordingState != MainVariables){ + if ((RecordingState != MainVariables) || (config->GetTime_Domain())){ MainRecording(); @@ -188,7 +188,7 @@ void CDiscAdjSinglezoneDriver::Run() { iteration->Iterate(output_container[ZONE_0], integration_container, geometry_container, solver_container, numerics_container, config_container, - surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); + surface_movement, grid_movement, FFDBox, ZONE_0, INST_0, false); /*--- Monitor the pseudo-time ---*/ @@ -313,7 +313,7 @@ void CDiscAdjSinglezoneDriver::SetRecording(unsigned short kind_recording){ void CDiscAdjSinglezoneDriver::SetAdj_ObjFunction(){ - bool time_stepping = config->GetTime_Marching() != STEADY; + bool time_stepping = config->GetTime_Domain(); unsigned long IterAvg_Obj = config->GetIter_Avg_Objective(); su2double seeding = 1.0; @@ -464,7 +464,7 @@ void CDiscAdjSinglezoneDriver::DirectRun(unsigned short kind_recording){ /*--- Iterate the direct solver ---*/ - direct_iteration->Iterate(direct_output, integration_container, geometry_container, solver_container, numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); + direct_iteration->Iterate(direct_output, integration_container, geometry_container, solver_container, numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, INST_0, false); /*--- Postprocess the direct solver ---*/ diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 841a0d6cb822..c284a00c1e8e 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -1288,7 +1288,10 @@ void CDriver::Solver_Restart(CSolver ***solver, CGeometry **geometry, solver[MESH_0][RAD_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); } if (fem) { - if (time_domain) val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; + if (time_domain) { + if (config->GetRestart()) val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; + else val_iter = SU2_TYPE::Int(config->GetUnst_AdjointIter())-1; + } solver[MESH_0][FEA_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); } if (fem_euler || fem_ns) { @@ -1333,7 +1336,8 @@ void CDriver::Solver_Restart(CSolver ***solver, CGeometry **geometry, if ((restart || restart_flow) && config->GetDeform_Mesh() && update_geo){ /*--- Always restart with the last state ---*/ - val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; + if (config->GetRestart()) val_iter = SU2_TYPE::Int(config->GetRestart_Iter())-1; + else val_iter = SU2_TYPE::Int(config->GetUnst_AdjointIter())-1; solver[MESH_0][MESH_SOL]->LoadRestart(geometry, solver, config, val_iter, update_geo); } @@ -2558,7 +2562,8 @@ void CDriver::Interface_Preprocessing(CConfig **config, CSolver***** solver, CGe "Use DEFORM_MESH=YES, and setup MARKER_DEFORM_MESH=(...)", CURRENT_FUNCTION); } interface_type = BOUNDARY_DISPLACEMENTS; - interface[donor][target] = new CDisplacementsInterface(nDim, 0); + if (!config[donor]->GetTime_Domain()) interface[donor][target] = new CDisplacementsInterface(nDim, 0); + else interface[donor][target] = new CDisplacementsInterface(2*nDim, 0); if (rank == MASTER_NODE) cout << "boundary displacements from the structural solver." << endl; } else if (fluid_donor && fluid_target) { @@ -3040,7 +3045,7 @@ void CFluidDriver::Run() { for (iZone = 0; iZone < nZone; iZone++) { config_container[iZone]->SetInnerIter(IntIter); - iteration_container[iZone][INST_0]->Iterate(output_container[iZone], integration_container, geometry_container, solver_container, numerics_container, config_container, surface_movement, grid_movement, FFDBox, iZone, INST_0); + iteration_container[iZone][INST_0]->Iterate(output_container[iZone], integration_container, geometry_container, solver_container, numerics_container, config_container, surface_movement, grid_movement, FFDBox, iZone, INST_0, false); } /*--- Check convergence in each zone --*/ @@ -3215,7 +3220,7 @@ void CTurbomachineryDriver::Run() { for (iZone = 0; iZone < nZone; iZone++) { iteration_container[iZone][INST_0]->Iterate(output_container[iZone], integration_container, geometry_container, solver_container, numerics_container, config_container, - surface_movement, grid_movement, FFDBox, iZone, INST_0); + surface_movement, grid_movement, FFDBox, iZone, INST_0, false); } for (iZone = 0; iZone < nZone; iZone++) { @@ -3458,7 +3463,7 @@ void CHBDriver::Run() { for (iInst = 0; iInst < nInstHB; iInst++) iteration_container[ZONE_0][iInst]->Iterate(output_container[ZONE_0], integration_container, geometry_container, solver_container, numerics_container, config_container, - surface_movement, grid_movement, FFDBox, ZONE_0, iInst); + surface_movement, grid_movement, FFDBox, ZONE_0, iInst, false); /*--- Update the convergence history file (serial and parallel computations). ---*/ diff --git a/SU2_CFD/src/drivers/CMultizoneDriver.cpp b/SU2_CFD/src/drivers/CMultizoneDriver.cpp index 5ab0163d21d3..8004bc539342 100644 --- a/SU2_CFD/src/drivers/CMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CMultizoneDriver.cpp @@ -307,7 +307,8 @@ void CMultizoneDriver::Run_GaussSeidel() { for (iZone = 0; iZone < nZone; iZone++){ /*--- In principle, the mesh does not need to be updated ---*/ - UpdateMesh = 0; + /*--- Default updated if config for zone specifies mesh deformation ---*/ + UpdateMesh = config_container[iZone]->GetDeform_Mesh(); /*--- Set the OuterIter ---*/ config_container[iZone]->SetOuterIter(iOuter_Iter); @@ -454,16 +455,23 @@ void CMultizoneDriver::Update() { for (auto jZone = 0u; jZone < nZone; jZone++){ /*--- The target zone is iZone ---*/ if (jZone != iZone){ + // Run update on structural solution to update relaxation on grid velocities before transferring data + if (config_container[jZone]->GetStructuralProblem()) { + iteration_container[jZone][INST_0]->Update(output_container[jZone], integration_container, geometry_container, + solver_container, numerics_container, config_container, + surface_movement, grid_movement, FFDBox, jZone, INST_0); + } UpdateMesh += Transfer_Data(jZone, iZone); } } /*--- If a mesh update is required due to the transfer of data ---*/ if (UpdateMesh > 0) DynamicMeshUpdate(iZone, TimeIter); - iteration_container[iZone][INST_0]->Update(output_container[iZone], integration_container, geometry_container, - solver_container, numerics_container, config_container, - surface_movement, grid_movement, FFDBox, iZone, INST_0); - + if (!config_container[iZone]->GetStructuralProblem()) { + iteration_container[iZone][INST_0]->Update(output_container[iZone], integration_container, geometry_container, + solver_container, numerics_container, config_container, + surface_movement, grid_movement, FFDBox, iZone, INST_0); + } /*--- Set the Convergence_FSI boolean to false for the next time step ---*/ for (unsigned short iSol = 0; iSol < MAX_SOLS-1; iSol++){ if (integration_container[iZone][INST_0][iSol] != nullptr){ diff --git a/SU2_CFD/src/interfaces/fsi/CDisplacementsInterface.cpp b/SU2_CFD/src/interfaces/fsi/CDisplacementsInterface.cpp index 4ea4f71c9e0b..7d255f95683d 100644 --- a/SU2_CFD/src/interfaces/fsi/CDisplacementsInterface.cpp +++ b/SU2_CFD/src/interfaces/fsi/CDisplacementsInterface.cpp @@ -43,12 +43,28 @@ void CDisplacementsInterface::GetDonor_Variable(CSolver *struct_solution, CGeome for (auto iVar = 0u; iVar < nVar; iVar++) Donor_Variable[iVar] = DisplacementDonor[iVar]; + + if (struct_config->GetTime_Domain()) { + auto VelocityDonor = struct_solution->GetNodes()->GetSolution_Vel_Pred(Point_Struct); + for (auto iVar = nVar/2; iVar < nVar; iVar++)//Assuming dynamic interface always has nVar = 2*nDim, 2D: 4, 3D: 6 + Donor_Variable[iVar] = VelocityDonor[iVar-nVar/2]; + } } void CDisplacementsInterface::SetTarget_Variable(CSolver *mesh_solver, CGeometry *flow_geometry, const CConfig *flow_config, unsigned long Marker_Flow, unsigned long Vertex_Flow, unsigned long Point_Mesh) { - /*--- Impose the boundary displacements ---*/ - mesh_solver->GetNodes()->SetBound_Disp(Point_Mesh,Target_Variable); + if (!flow_config->GetTime_Domain()) { + /*--- Impose the boundary displacements ---*/ + mesh_solver->GetNodes()->SetBound_Disp(Point_Mesh,Target_Variable); + } else { + /*--- Impose the boundary displacements ---*/ + for (auto iVar = 0u; iVar < nVar/2; iVar++) + mesh_solver->GetNodes()->SetBound_Disp(Point_Mesh,iVar,Target_Variable[iVar]); + + /*--- Impose the boundary velocities ---*/ + for (auto iVar = nVar/2; iVar < nVar; iVar++) + mesh_solver->GetNodes()->SetBound_Vel(Point_Mesh,iVar-nVar/2,Target_Variable[iVar]); + } } diff --git a/SU2_CFD/src/iteration/CAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CAdjFluidIteration.cpp index 0484d9d7654d..e24e447fc71a 100644 --- a/SU2_CFD/src/iteration/CAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CAdjFluidIteration.cpp @@ -133,7 +133,7 @@ void CAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integratio void CAdjFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, - CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) { + CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst, bool CrossTerm) { switch (config[val_iZone]->GetKind_Solver()) { case ADJ_EULER: config[val_iZone]->SetGlobalParam(ADJ_EULER, RUNTIME_ADJFLOW_SYS); diff --git a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp index d889f036b4ce..2355451949e9 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp @@ -161,7 +161,7 @@ void CDiscAdjFEAIteration::LoadDynamic_Solution(CGeometry**** geometry, CSolver* void CDiscAdjFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** volume_grid_movement, - CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) { + CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst, bool CrossTerm) { bool dynamic = (config[val_iZone]->GetTime_Domain()); /*--- Extract the adjoints of the conservative input variables and store them for the next iteration ---*/ @@ -207,7 +207,7 @@ void CDiscAdjFEAIteration::SetRecording(COutput* output, CIntegration**** integr /*--- Run one iteration while tape is passive - this clears all indices ---*/ fem_iteration->Iterate(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, - FFDBox, val_iZone, val_iInst); + FFDBox, val_iZone, val_iInst, false); } /*--- Prepare for recording ---*/ @@ -236,7 +236,7 @@ void CDiscAdjFEAIteration::SetRecording(COutput* output, CIntegration**** integr /*--- Run the direct iteration ---*/ fem_iteration->Iterate(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, - FFDBox, val_iZone, val_iInst); + FFDBox, val_iZone, val_iInst, false); config[val_iZone]->SetTimeIter(TimeIter); @@ -280,11 +280,10 @@ void CDiscAdjFEAIteration::RegisterInput(CSolver***** solver, CGeometry**** geom /*--- Register topology optimization densities (note direct solver) ---*/ solver[iZone][iInst][MESH_0][FEA_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]); - - /*--- Register mesh coordinates for geometric sensitivities ---*/ - - geometry[iZone][iInst][MESH_0]->RegisterCoordinates(config[iZone]); } + /*--- Register mesh coordinates for geometric sensitivities ---*/ + + geometry[iZone][iInst][MESH_0]->RegisterCoordinates(config[iZone]); } void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** geometry, CNumerics****** numerics, @@ -300,6 +299,7 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** ge bool nonlinear = config[iZone]->GetGeometricConditions() == LARGE_DEFORMATIONS; bool de_effects = config[iZone]->GetDE_Effects() && nonlinear; bool element_based = dir_solver->IsElementBased() && nonlinear; + bool dynamic = config[iZone]->GetTime_Domain(); for (unsigned short iProp = 0; iProp < config[iZone]->GetnElasticityMod(); iProp++) { su2double E = adj_solver->GetVal_Young(iProp); @@ -368,6 +368,7 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** ge if (fsi) { /*--- Set relation between solution and predicted displacements, which are the transferred ones. ---*/ dir_solver->PredictStruct_Displacement(structural_geometry, config[iZone]); + if (dynamic) dir_solver->PredictStruct_Velocity(structural_geometry, config[iZone]); } /*--- MPI dependencies. ---*/ @@ -375,10 +376,8 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** ge dir_solver->InitiateComms(structural_geometry, config[iZone], SOLUTION_FEA); dir_solver->CompleteComms(structural_geometry, config[iZone], SOLUTION_FEA); - if (kind_recording == MESH_COORDS) { - structural_geometry->InitiateComms(structural_geometry, config[iZone], COORDINATES); - structural_geometry->CompleteComms(structural_geometry, config[iZone], COORDINATES); - } + structural_geometry->InitiateComms(structural_geometry, config[iZone], COORDINATES); + structural_geometry->CompleteComms(structural_geometry, config[iZone], COORDINATES); /*--- Topology optimization dependencies. ---*/ diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index 22c2b384c5de..b34c48dd6b45 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -87,6 +87,10 @@ void CDiscAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integr geometry[val_iZone][val_iInst][iMesh]->nodes->SetCoord_n(); geometry[val_iZone][val_iInst][iMesh]->nodes->SetCoord_n1(); } + if (config[val_iZone]->GetDynamic_Grid()) { + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume_n(); + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume_nM1(); + } } } if (dual_time) { @@ -106,6 +110,9 @@ void CDiscAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integr if (grid_IsMoving) { geometry[val_iZone][val_iInst][iMesh]->nodes->SetCoord_n(); } + if (config[val_iZone]->GetDynamic_Grid()) { + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume_n(); + } } } @@ -126,6 +133,16 @@ void CDiscAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integr Afterwards the GridVelocity is computed based on the Coordinates. ---*/ + /*--- Temporarily store the loaded volumes in to old containers ---*/ + if (config[val_iZone]->GetDynamic_Grid()) { + for (iMesh=0; iMesh<=config[val_iZone]->GetnMGLevels();iMesh++) { + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume_Old(); + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume_n_Old(); + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume_nM1_Old(); + } + } + + /*-- Load mesh solver ---*/ if (config[val_iZone]->GetDeform_Mesh()) { solver[val_iZone][val_iInst][MESH_0][MESH_SOL]->LoadRestart( geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone], Direct_Iter, true); @@ -136,6 +153,26 @@ void CDiscAdjFluidIteration::Preprocess(COutput* output, CIntegration**** integr LoadUnsteady_Solution(geometry, solver, config, val_iInst, val_iZone, Direct_Iter - 1); } else { LoadUnsteady_Solution(geometry, solver, config, val_iInst, val_iZone, Direct_Iter - 2); + + /*--- Set volumes into correct containers ---*/ + if (config[val_iZone]->GetDynamic_Grid()) { + for (iMesh=0; iMesh<=config[val_iZone]->GetnMGLevels();iMesh++) { + /*--- If negative iteration number, set default ---*/ + if (Direct_Iter - 2 < 0) { + for(iPoint=0; iPointGetnPoint();iPoint++) { + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume(iPoint,0.0); + } + } + + /*--- Set currently loaded volume to Volume_nM1 ---*/ + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume_n(); + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume_nM1(); + + /*--- Set Volume_n and Volume from old containers ---*/ + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume_n_from_OldnM1(); + geometry[val_iZone][val_iInst][iMesh]->nodes->SetVolume_from_Oldn(); + } + } } /*--- Temporarily store the loaded solution in the Solution_Old array ---*/ @@ -342,25 +379,32 @@ void CDiscAdjFluidIteration::LoadUnsteady_Solution(CGeometry**** geometry, CSolv void CDiscAdjFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** volume_grid_movement, - CFreeFormDefBox*** FFDBox, unsigned short iZone, unsigned short iInst) { + CFreeFormDefBox*** FFDBox, unsigned short iZone, unsigned short iInst, bool CrossTerm) { bool frozen_visc = config[iZone]->GetFrozen_Visc_Disc(); bool heat = config[iZone]->GetWeakly_Coupled_Heat(); /*--- Extract the adjoints of the conservative input variables and store them for the next iteration ---*/ if (config[iZone]->GetFluidProblem()) { - solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone]); + solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm); solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]); + + /*--- If mesh deformation defined in config, extract here as well ---*/ + if (config[iZone]->GetDeform_Mesh()) { + solver[iZone][iInst][MESH_0][ADJMESH_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone]); + + solver[iZone][iInst][MESH_0][ADJMESH_SOL]->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]); + } } if (turbulent && !frozen_visc) { - solver[iZone][iInst][MESH_0][ADJTURB_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone]); + solver[iZone][iInst][MESH_0][ADJTURB_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm); } if (heat) { - solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone]); + solver[iZone][iInst][MESH_0][ADJHEAT_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm); } if (config[iZone]->AddRadiation()) { - solver[iZone][iInst][MESH_0][ADJRAD_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone]); + solver[iZone][iInst][MESH_0][ADJRAD_SOL]->ExtractAdjoint_Solution(geometry[iZone][iInst][MESH_0], config[iZone], CrossTerm); solver[iZone][iInst][MESH_0][ADJRAD_SOL]->ExtractAdjoint_Variables(geometry[iZone][iInst][MESH_0], config[iZone]); } @@ -407,6 +451,16 @@ void CDiscAdjFluidIteration::RegisterInput(CSolver***** solver, CGeometry**** ge solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]); solver[iZone][iInst][MESH_0][ADJFLOW_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]); + + /*--- If mesh deformation defined in config, register here as well ---*/ + if (config[iZone]->GetDeform_Mesh()) { + + /*--- Undeformed mesh coordinates ---*/ + solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterSolution(geometry[iZone][iInst][MESH_0], config[iZone]); + + /*--- Boundary displacements ---*/ + solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterVariables(geometry[iZone][iInst][MESH_0], config[iZone]); + } } if (turbulent && !frozen_visc) { diff --git a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp index 041ea3fcc06a..1631259b7bac 100644 --- a/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjHeatIteration.cpp @@ -169,7 +169,7 @@ void CDiscAdjHeatIteration::LoadUnsteady_Solution(CGeometry**** geometry, CSolve void CDiscAdjHeatIteration::Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** volume_grid_movement, - CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) { + CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst, bool CrossTerm) { solver[val_iZone][val_iInst][MESH_0][ADJHEAT_SOL]->ExtractAdjoint_Solution(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); } diff --git a/SU2_CFD/src/iteration/CFEAIteration.cpp b/SU2_CFD/src/iteration/CFEAIteration.cpp index 1fd3f093a60f..f59c5dcb16cd 100644 --- a/SU2_CFD/src/iteration/CFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CFEAIteration.cpp @@ -31,7 +31,7 @@ void CFEAIteration::Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst) { + unsigned short val_iInst, bool CrossTerm) { bool StopCalc = false; unsigned long IntIter = 0; @@ -206,7 +206,6 @@ void CFEAIteration::Update(COutput* output, CIntegration**** integration, CGeome const su2double Physical_t = (TimeIter + 1) * Physical_dt; if (Physical_t >= config[val_iZone]->GetTotal_DynTime()) integration[val_iZone][val_iInst][FEA_SOL]->SetConvergence(true); - } else if (fsi) { /*--- For FSI problems, output the relaxed result, which is the one transferred into the fluid domain (for restart * purposes) ---*/ @@ -223,6 +222,7 @@ void CFEAIteration::Predictor(COutput* output, CIntegration**** integration, CGe CSolver* feaSolver = solver[val_iZone][val_iInst][MESH_0][FEA_SOL]; feaSolver->PredictStruct_Displacement(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); + if (config[val_iZone]->GetTime_Domain()) feaSolver->PredictStruct_Velocity(geometry[val_iZone][val_iInst][MESH_0], config[val_iZone]); } void CFEAIteration::Relaxation(COutput* output, CIntegration**** integration, CGeometry**** geometry, @@ -266,7 +266,7 @@ void CFEAIteration::Solve(COutput* output, CIntegration**** integration, CGeomet unsigned short val_iInst) { /*------------------ Structural subiteration ----------------------*/ Iterate(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox, val_iZone, - val_iInst); + val_iInst, false); if (multizone && !config[val_iZone]->GetTime_Domain()) { Output(output, geometry, solver, config, config[val_iZone]->GetOuterIter(), false, val_iZone, val_iInst); diff --git a/SU2_CFD/src/iteration/CFEMFluidIteration.cpp b/SU2_CFD/src/iteration/CFEMFluidIteration.cpp index cabea1e48c82..4e2049fd4420 100644 --- a/SU2_CFD/src/iteration/CFEMFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFEMFluidIteration.cpp @@ -44,7 +44,7 @@ void CFEMFluidIteration::Preprocess(COutput* output, CIntegration**** integratio void CFEMFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, - CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) { + CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst, bool CrossTerm) { /*--- Update global parameters ---*/ if (config[val_iZone]->GetKind_Solver() == FEM_EULER || config[val_iZone]->GetKind_Solver() == DISC_ADJ_FEM_EULER) diff --git a/SU2_CFD/src/iteration/CFluidIteration.cpp b/SU2_CFD/src/iteration/CFluidIteration.cpp index b53b35f7bc32..d4a5078df12c 100644 --- a/SU2_CFD/src/iteration/CFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CFluidIteration.cpp @@ -36,11 +36,12 @@ void CFluidIteration::Preprocess(COutput* output, CIntegration**** integration, bool fsi = config[val_iZone]->GetFSI_Simulation(); unsigned long OuterIter = config[val_iZone]->GetOuterIter(); + bool adjoint = (config[val_iZone]->GetContinuous_Adjoint() || config[val_iZone]->GetDiscrete_Adjoint()); /*--- Set the initial condition for FSI problems with subiterations ---*/ /*--- This is done only in the first block subiteration.---*/ /*--- From then on, the solver reuses the partially converged solution obtained in the previous subiteration ---*/ - if (fsi && (OuterIter == 0)) { + if (fsi && !adjoint && (OuterIter == 0)) { solver[val_iZone][val_iInst][MESH_0][FLOW_SOL]->SetInitialCondition( geometry[val_iZone][val_iInst], solver[val_iZone][val_iInst], config[val_iZone], TimeIter); } @@ -55,7 +56,7 @@ void CFluidIteration::Preprocess(COutput* output, CIntegration**** integration, void CFluidIteration::Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, - CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst) { + CFreeFormDefBox*** FFDBox, unsigned short val_iZone, unsigned short val_iInst, bool CrossTerm) { unsigned long InnerIter, TimeIter; const bool unsteady = (config[val_iZone]->GetTime_Marching() == DT_STEPPING_1ST) || @@ -286,7 +287,7 @@ void CFluidIteration::Solve(COutput* output, CIntegration**** integration, CGeom /*--- Run a single iteration of the solver ---*/ Iterate(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox, val_iZone, - INST_0); + INST_0, false); /*--- Monitor the pseudo-time ---*/ StopCalc = Monitor(output, integration, geometry, solver, numerics, config, surface_movement, grid_movement, FFDBox, diff --git a/SU2_CFD/src/iteration/CHeatIteration.cpp b/SU2_CFD/src/iteration/CHeatIteration.cpp index 50fcd31140b5..2396b0dbb4bd 100644 --- a/SU2_CFD/src/iteration/CHeatIteration.cpp +++ b/SU2_CFD/src/iteration/CHeatIteration.cpp @@ -31,7 +31,7 @@ void CHeatIteration::Iterate(COutput* output, CIntegration**** integration, CGeometry**** geometry, CSolver***** solver, CNumerics****** numerics, CConfig** config, CSurfaceMovement** surface_movement, CVolumetricMovement*** grid_movement, CFreeFormDefBox*** FFDBox, unsigned short val_iZone, - unsigned short val_iInst) { + unsigned short val_iInst, bool CrossTerm) { /*--- Update global parameters ---*/ config[val_iZone]->SetGlobalParam(HEAT_EQUATION, RUNTIME_HEAT_SYS); diff --git a/SU2_CFD/src/output/CElasticityOutput.cpp b/SU2_CFD/src/output/CElasticityOutput.cpp index b5a8ae3bd243..bbff044e138e 100644 --- a/SU2_CFD/src/output/CElasticityOutput.cpp +++ b/SU2_CFD/src/output/CElasticityOutput.cpp @@ -72,6 +72,10 @@ CElasticityOutput::CElasticityOutput(CConfig *config, unsigned short nDim) : COu requestedVolumeFields.emplace_back("COORDINATES"); requestedVolumeFields.emplace_back("SOLUTION"); requestedVolumeFields.emplace_back("STRESS"); + if (dynamic) { + requestedVolumeFields.emplace_back("VELOCITY"); + requestedVolumeFields.emplace_back("ACCELERATION"); + } if (config->GetTopology_Optimization()) requestedVolumeFields.emplace_back("TOPOLOGY"); nRequestedVolumeFields = requestedVolumeFields.size(); } diff --git a/SU2_CFD/src/output/CFlowCompOutput.cpp b/SU2_CFD/src/output/CFlowCompOutput.cpp index e59202bf7aef..98dbef3a7687 100644 --- a/SU2_CFD/src/output/CFlowCompOutput.cpp +++ b/SU2_CFD/src/output/CFlowCompOutput.cpp @@ -35,7 +35,7 @@ CFlowCompOutput::CFlowCompOutput(CConfig *config, unsigned short nDim) : CFlowOu turb_model = config->GetKind_Turb_Model(); lastInnerIter = curInnerIter; - gridMovement = config->GetGrid_Movement(); + gridMovement = (config->GetGrid_Movement() || config->GetDynamic_Grid()); /*--- Set the default history fields if nothing is set in the config file ---*/ @@ -58,7 +58,7 @@ CFlowCompOutput::CFlowCompOutput(CConfig *config, unsigned short nDim) : CFlowOu requestedVolumeFields.emplace_back("COORDINATES"); requestedVolumeFields.emplace_back("SOLUTION"); requestedVolumeFields.emplace_back("PRIMITIVE"); - if (config->GetGrid_Movement()) requestedVolumeFields.emplace_back("GRID_VELOCITY"); + if (gridMovement) requestedVolumeFields.emplace_back("GRID_VELOCITY"); nRequestedVolumeFields = requestedVolumeFields.size(); } @@ -313,7 +313,7 @@ void CFlowCompOutput::SetVolumeOutputFields(CConfig *config){ } // Grid velocity - if (config->GetGrid_Movement()){ + if (gridMovement) { AddVolumeOutput("GRID_VELOCITY-X", "Grid_Velocity_x", "GRID_VELOCITY", "x-component of the grid velocity vector"); AddVolumeOutput("GRID_VELOCITY-Y", "Grid_Velocity_y", "GRID_VELOCITY", "y-component of the grid velocity vector"); if (nDim == 3 ) @@ -467,7 +467,7 @@ void CFlowCompOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolv break; } - if (config->GetGrid_Movement()){ + if (gridMovement){ SetVolumeOutputValue("GRID_VELOCITY-X", iPoint, Node_Geo->GetGridVel(iPoint)[0]); SetVolumeOutputValue("GRID_VELOCITY-Y", iPoint, Node_Geo->GetGridVel(iPoint)[1]); if (nDim == 3) diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp index 927cb743c8ff..518343381b78 100644 --- a/SU2_CFD/src/output/CFlowIncOutput.cpp +++ b/SU2_CFD/src/output/CFlowIncOutput.cpp @@ -38,7 +38,7 @@ CFlowIncOutput::CFlowIncOutput(CConfig *config, unsigned short nDim) : CFlowOutp heat = config->GetEnergy_Equation(); weakly_coupled_heat = config->GetWeakly_Coupled_Heat(); - + gridMovement = (config->GetGrid_Movement() || config->GetDynamic_Grid()); /*--- Set the default history fields if nothing is set in the config file ---*/ if (nRequestedHistoryFields == 0){ @@ -60,7 +60,7 @@ CFlowIncOutput::CFlowIncOutput(CConfig *config, unsigned short nDim) : CFlowOutp requestedVolumeFields.emplace_back("COORDINATES"); requestedVolumeFields.emplace_back("SOLUTION"); requestedVolumeFields.emplace_back("PRIMITIVE"); - if (config->GetGrid_Movement()) requestedVolumeFields.emplace_back("GRID_VELOCITY"); + if (gridMovement) requestedVolumeFields.emplace_back("GRID_VELOCITY"); nRequestedVolumeFields = requestedVolumeFields.size(); } @@ -383,7 +383,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){ AddVolumeOutput("P1-RAD", "Radiative_Energy(P1)", "SOLUTION", "Radiative Energy"); // Grid velocity - if (config->GetGrid_Movement()){ + if (gridMovement){ AddVolumeOutput("GRID_VELOCITY-X", "Grid_Velocity_x", "GRID_VELOCITY", "x-component of the grid velocity vector"); AddVolumeOutput("GRID_VELOCITY-Y", "Grid_Velocity_y", "GRID_VELOCITY", "y-component of the grid velocity vector"); if (nDim == 3 ) diff --git a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp index 1621d22f12f0..a4e9f8862e8e 100644 --- a/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp @@ -350,15 +350,13 @@ void CDiscAdjFEASolver::RegisterSolution(CGeometry *geometry, CConfig *config){ /*--- Register acceleration (u'') and velocity (u') at time step n ---*/ - direct_solver->GetNodes()->RegisterSolution_Accel(input); - direct_solver->GetNodes()->RegisterSolution_Vel(input); + direct_solver->GetNodes()->RegisterSolution_Accel(input, push_index); + direct_solver->GetNodes()->RegisterSolution_Vel(input, push_index); /*--- Register solution (u), acceleration (u'') and velocity (u') at time step n-1 ---*/ - - direct_solver->GetNodes()->Register_femSolution_time_n(); - direct_solver->GetNodes()->RegisterSolution_Accel_time_n(); - direct_solver->GetNodes()->RegisterSolution_Vel_time_n(); - + direct_solver->GetNodes()->Register_femSolution_time_n(input, push_index); + direct_solver->GetNodes()->RegisterSolution_Accel_time_n(input, push_index); + direct_solver->GetNodes()->RegisterSolution_Vel_time_n(input, push_index); } } @@ -456,8 +454,8 @@ void CDiscAdjFEASolver::RegisterOutput(CGeometry *geometry, CConfig *config){ if (dynamic) { /*--- Register acceleration (u'') and velocity (u') at time step n ---*/ - direct_solver->GetNodes()->RegisterSolution_Accel(input); - direct_solver->GetNodes()->RegisterSolution_Vel(input); + direct_solver->GetNodes()->RegisterSolution_Accel(input, push_index); + direct_solver->GetNodes()->RegisterSolution_Vel(input, push_index); } } @@ -549,7 +547,7 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co /*--- Extract the adjoint solution ---*/ - if(config->GetMultizone_Problem()) { + if(multizone) { direct_solver->GetNodes()->GetAdjointSolution_LocalIndex(iPoint,Solution); } else { @@ -575,7 +573,12 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co /*--- Extract the adjoint acceleration solution u'' ---*/ + if(multizone) { + direct_solver->GetNodes()->GetAdjointSolution_Accel_LocalIndex(iPoint,Solution_Accel); + } + else { direct_solver->GetNodes()->GetAdjointSolution_Accel(iPoint,Solution_Accel); + } /*--- Store the adjoint acceleration solution u'' ---*/ @@ -591,8 +594,12 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co for (iPoint = 0; iPoint < nPoint; iPoint++){ /*--- Extract the adjoint velocity solution u'' ---*/ - + if(multizone) { + direct_solver->GetNodes()->GetAdjointSolution_Vel_LocalIndex(iPoint,Solution_Vel); + } + else { direct_solver->GetNodes()->GetAdjointSolution_Vel(iPoint,Solution_Vel); + } /*--- Store the adjoint velocity solution u'' ---*/ @@ -604,8 +611,12 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co for (iPoint = 0; iPoint < nPoint; iPoint++){ /*--- Extract the adjoint solution at time n ---*/ - + if(multizone) { + direct_solver->GetNodes()->GetAdjointSolution_time_n_LocalIndex(iPoint,Solution); + } + else { direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); + } /*--- Store the adjoint solution at time n ---*/ @@ -616,8 +627,12 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co for (iPoint = 0; iPoint < nPoint; iPoint++){ /*--- Extract the adjoint acceleration solution u'' at time n ---*/ - + if(multizone) { + direct_solver->GetNodes()->GetAdjointSolution_Accel_time_n_LocalIndex(iPoint,Solution_Accel); + } + else { direct_solver->GetNodes()->GetAdjointSolution_Accel_time_n(iPoint,Solution_Accel); + } /*--- Store the adjoint acceleration solution u'' at time n---*/ @@ -629,8 +644,12 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co for (iPoint = 0; iPoint < nPoint; iPoint++){ /*--- Extract the adjoint velocity solution u' at time n ---*/ - + if(multizone) { + direct_solver->GetNodes()->GetAdjointSolution_Vel_time_n_LocalIndex(iPoint,Solution_Vel); + } + else { direct_solver->GetNodes()->GetAdjointSolution_Vel_time_n(iPoint,Solution_Vel); + } /*--- Store the adjoint velocity solution u' at time n ---*/ @@ -679,8 +698,7 @@ void CDiscAdjFEASolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig *c unsigned short iVar; bool local_index = config->GetMultizone_Problem(); - /*--- Extract the adjoint values of the farfield values ---*/ - + /*--- Extract the adjoint values of the material properties ---*/ if (KindDirect_Solver == RUNTIME_FEA_SYS){ if (local_index) { @@ -762,6 +780,9 @@ void CDiscAdjFEASolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config){ } for (iVar = 0; iVar < nVar; iVar++){ Solution_Vel[iVar] = nodes->GetSolution_Vel(iPoint,iVar); + if (deform_mesh){ + Solution_Vel[iVar] += nodes->GetSourceTerm_VelAdjoint(iPoint,iVar); + } } for (iVar = 0; iVar < nVar; iVar++){ Solution[iVar] += nodes->GetDynamic_Derivative_n(iPoint,iVar); @@ -807,6 +828,8 @@ void CDiscAdjFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_cont void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*){ + bool time_domain = config->GetTime_Domain(); + bool multizone = config->GetMultizone_Problem(); unsigned short iVar; for (iVar = 0; iVar < nMPROP; iVar++){ @@ -840,7 +863,7 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo su2double Sensitivity; - if(config->GetMultizone_Problem()) { + if(multizone) { Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim); } else { @@ -848,11 +871,75 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo /*--- Set the index manually to zero. ---*/ AD::ResetInput(Coord[iDim]); } - - nodes->SetSensitivity(iPoint, iDim, Sensitivity); + if (!time_domain || multizone) { + nodes->SetSensitivity(iPoint, iDim, Sensitivity); + } else { + nodes->SetSensitivity(iPoint, iDim, nodes->GetSensitivity(iPoint, iDim) + Sensitivity); + } } } SetSurface_Sensitivity(geometry, config); + + // Temporary Output + if (config->GetAdvanced_FEAElementBased()) { + if (rank == MASTER_NODE) { + unsigned short iVar; + + /*--- Header of the temporary adjoint output file ---*/ + ofstream myfile_res; + myfile_res.open("Results_Reverse_Adjoint.txt", ios::app); + myfile_res.precision(15); + myfile_res << config->GetTimeIter() << "\t"; + switch (config->GetKind_ObjFunc()) { + case REFERENCE_NODE: + myfile_res << scientific << direct_solver->GetTotal_OFRefNode() << "\t"; + break; + default: + myfile_res << scientific << 0.0 << "\t"; + break; + } + for (iVar = 0; iVar < nMPROP; iVar++) + myfile_res << scientific << Total_Sens_E[iVar] << "\t"; + for (iVar = 0; iVar < nDV; iVar++) { + myfile_res << scientific << Total_Sens_DV[iVar] << "\t"; + } + myfile_res << endl; + myfile_res.close(); + + + /*--- Header of the temporary dv sensitivity output file ---*/ + ofstream myfile2_res; + bool outputDVFEA = false; + + switch (config->GetDV_FEA()) { + case YOUNG_MODULUS: + myfile2_res.open("grad_young.opt"); + outputDVFEA = true; + break; + default: + outputDVFEA = false; + break; + } + + if (outputDVFEA) { + unsigned short iDV; + + myfile2_res << "INDEX" + << "\t" + << "GRAD" << endl; + + myfile2_res.precision(15); + + for (iDV = 0; iDV < nDV; iDV++) { + myfile2_res << iDV; + myfile2_res << "\t"; + myfile2_res << scientific << Total_Sens_DV[iDV]; + myfile2_res << endl; + } + } + myfile2_res.close(); + } + } } void CDiscAdjFEASolver::SetSurface_Sensitivity(CGeometry *geometry, CConfig *config){ diff --git a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp index de7eebe5d27e..a64dedfc973c 100644 --- a/SU2_CFD/src/solvers/CDiscAdjSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjSolver.cpp @@ -261,10 +261,10 @@ void CDiscAdjSolver::RegisterSolution(CGeometry *geometry, CConfig *config) { direct_solver->GetNodes()->RegisterSolution(input, push_index); if (time_n_needed) - direct_solver->GetNodes()->RegisterSolution_time_n(); + direct_solver->GetNodes()->RegisterSolution_time_n(push_index); if (time_n1_needed) - direct_solver->GetNodes()->RegisterSolution_time_n1(); + direct_solver->GetNodes()->RegisterSolution_time_n1(push_index); } void CDiscAdjSolver::RegisterVariables(CGeometry *geometry, CConfig *config, bool reset) { @@ -475,7 +475,7 @@ void CDiscAdjSolver::SetAdj_ObjFunc(CGeometry *geometry, CConfig *config) { } } -void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config){ +void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm){ const bool time_n1_needed = config->GetTime_Marching() == DT_STEPPING_2ND; const bool time_n_needed = (config->GetTime_Marching() == DT_STEPPING_1ST) || time_n1_needed; @@ -500,7 +500,7 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi /*--- Extract the adjoint solution ---*/ - if(config->GetMultizone_Problem()) { + if(multizone) { direct_solver->GetNodes()->GetAdjointSolution_LocalIndex(iPoint,Solution); } else { @@ -529,11 +529,16 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi /*--- Extract the adjoint solution at time n ---*/ + if(multizone) { + direct_solver->GetNodes()->GetAdjointSolution_time_n_LocalIndex(iPoint,Solution); + } + else { direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution); + } /*--- Store the adjoint solution at time n ---*/ - nodes->Set_Solution_time_n(iPoint,Solution); + if (!CrossTerm) nodes->Set_Solution_time_n(iPoint,Solution); } } @@ -542,11 +547,16 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi /*--- Extract the adjoint solution at time n-1 ---*/ + if(multizone) { + direct_solver->GetNodes()->GetAdjointSolution_time_n1_LocalIndex(iPoint,Solution); + } + else { direct_solver->GetNodes()->GetAdjointSolution_time_n1(iPoint,Solution); + } /*--- Store the adjoint solution at time n-1 ---*/ - nodes->Set_Solution_time_n1(iPoint,Solution); + if (!CrossTerm) nodes->Set_Solution_time_n1(iPoint,Solution); } } @@ -623,7 +633,7 @@ void CDiscAdjSolver::ExtractAdjoint_Geometry(CGeometry *geometry, CConfig *confi // unsigned short iVar; unsigned long iPoint; - + bool multizone = config->GetMultizone_Problem(); /*--- Set Residuals to zero ---*/ // for (iVar = 0; iVar < nVar; iVar++){ @@ -639,7 +649,7 @@ void CDiscAdjSolver::ExtractAdjoint_Geometry(CGeometry *geometry, CConfig *confi /*--- Extract the adjoint solution ---*/ - if (config->GetMultizone_Problem()) + if (multizone) geometry->nodes->GetAdjointCoord_LocalIndex(iPoint, Solution_Geometry); else geometry->nodes->GetAdjointCoord(iPoint, Solution_Geometry); @@ -693,7 +703,7 @@ void CDiscAdjSolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config) { bool dual_time = (config->GetTime_Marching() == DT_STEPPING_1ST || config->GetTime_Marching() == DT_STEPPING_2ND); - + bool multizone = config->GetMultizone_Problem(); unsigned short iVar; unsigned long iPoint; @@ -706,7 +716,7 @@ void CDiscAdjSolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config) { Solution[iVar] += nodes->GetDual_Time_Derivative(iPoint,iVar); } } - if(config->GetMultizone_Problem()) { + if(multizone) { direct_solver->GetNodes()->SetAdjointSolution_LocalIndex(iPoint,Solution); } else { @@ -745,7 +755,7 @@ void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolve unsigned long iPoint; unsigned short iDim; su2double *Coord, Sensitivity, eps; - + bool multizone = config->GetMultizone_Problem(); bool time_stepping = (config->GetTime_Marching() != STEADY); for (iPoint = 0; iPoint < nPoint; iPoint++) { @@ -753,7 +763,7 @@ void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolve for (iDim = 0; iDim < nDim; iDim++) { - if(config->GetMultizone_Problem()) { + if(multizone) { Sensitivity = geometry->nodes->GetAdjointSolution(iPoint, iDim); } else { @@ -771,7 +781,7 @@ void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolve if ( geometry->nodes->GetSharpEdge_Distance(iPoint) < config->GetAdjSharp_LimiterCoeff()*eps ) Sensitivity = 0.0; } - if (!time_stepping) { + if (!time_stepping || multizone) { nodes->SetSensitivity(iPoint,iDim, Sensitivity); } else { nodes->SetSensitivity(iPoint, iDim, nodes->GetSensitivity(iPoint,iDim) + Sensitivity); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index fe10390bf57c..c1ed2e95e545 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -2240,6 +2240,26 @@ void CFEASolver::BC_Deforming(CGeometry *geometry, CNumerics *numerics, const CC } +void CFEASolver::BC_Velocity(CGeometry *geometry, CNumerics *numerics, const CConfig *config, unsigned short val_marker){ + + for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) { + + /*--- Get node index ---*/ + auto iNode = geometry->vertex[val_marker][iVertex]->GetNode(); + + /*--- Retrieve the boundary velocity ---*/ + su2double Vel[MAXNVAR] = {0.0}; + for (unsigned short iDim = 0; iDim < nDim; iDim++) + Vel[iDim] = nodes->GetBound_Vel(iNode,iDim); + + /*--- Set and enforce solution ---*/ + LinSysSol.SetBlock(iNode, Vel); + Jacobian.EnforceSolutionAtNode(iNode, Vel, LinSysRes); + + } + +} + su2double CFEASolver::Compute_LoadCoefficient(su2double CurrentTime, su2double RampTime, const CConfig *config){ su2double LoadCoeff = 1.0; @@ -2290,14 +2310,14 @@ su2double CFEASolver::Compute_LoadCoefficient(su2double CurrentTime, su2double R LoadCoeff = min(LoadCoeff,1.0); } - else if (Sine_Load){ + if (Sine_Load){ //Allow sine load with initial ramp loading as well /*--- Retrieve amplitude, frequency (Hz) and phase (rad) ---*/ su2double SineAmp = config->GetLoad_Sine()[0]; su2double SineFreq = config->GetLoad_Sine()[1]; su2double SinePhase = config->GetLoad_Sine()[2]; - LoadCoeff = SineAmp * sin(2*PI_NUMBER*SineFreq*CurrentTime + SinePhase); + LoadCoeff = LoadCoeff * SineAmp * sin(2*PI_NUMBER*SineFreq*CurrentTime + SinePhase); } /*--- Add possibility to release the load after the ramp---*/ @@ -2781,6 +2801,19 @@ void CFEASolver::PredictStruct_Displacement(CGeometry *geometry, CConfig *config } +void CFEASolver::PredictStruct_Velocity(CGeometry *geometry, CConfig *config) { + + /*--- To nPointDomain: we need to communicate the predicted solution after setting it. ---*/ + SU2_OMP_PARALLEL_(for schedule(static,omp_chunk_size)) + for (unsigned long iPoint=0; iPoint < nPointDomain; iPoint++) { + nodes->SetSolution_Vel_Pred(iPoint); + } + + InitiateComms(geometry, config, SOLUTION_VEL_PRED); + CompleteComms(geometry, config, SOLUTION_VEL_PRED); + +} + void CFEASolver::ComputeAitken_Coefficient(CGeometry *geometry, CConfig *config, unsigned long iOuterIter) { unsigned long iPoint, iDim; @@ -2890,6 +2923,10 @@ void CFEASolver::SetAitken_Relaxation(CGeometry *geometry, CConfig *config) { /*--- Set calculated solution as the old solution (needed for dynamic Aitken relaxation) ---*/ nodes->SetSolution_Old(iPoint, dispCalc); + /*--- Set predicted velocity to update in multizone iterations ---*/ + + if (config->GetTime_Domain()) nodes->SetSolution_Vel_Pred(iPoint); + /*--- Apply the Aitken relaxation ---*/ for (unsigned short iDim=0; iDim < nDim; iDim++) { dispPred[iDim] = (1.0 - WAitken)*dispPred[iDim] + WAitken*dispCalc[iDim]; @@ -3294,15 +3331,16 @@ void CFEASolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *c for (unsigned short iVar = 0; iVar < nVar; iVar++) { nodes->SetSolution(iPoint_Local, iVar, Sol[iVar]); if (dynamic) { - nodes->Set_Solution_time_n(iPoint_Local, iVar, Sol[iVar]); + if (!discrete_adjoint) nodes->Set_Solution_time_n(iPoint_Local, iVar, Sol[iVar]); nodes->SetSolution_Vel(iPoint_Local, iVar, Sol[iVar+nVar]); - nodes->SetSolution_Vel_time_n(iPoint_Local, iVar, Sol[iVar+nVar]); + if (!discrete_adjoint) nodes->SetSolution_Vel_time_n(iPoint_Local, iVar, Sol[iVar+nVar]); nodes->SetSolution_Accel(iPoint_Local, iVar, Sol[iVar+2*nVar]); - nodes->SetSolution_Accel_time_n(iPoint_Local, iVar, Sol[iVar+2*nVar]); + if (!discrete_adjoint) nodes->SetSolution_Accel_time_n(iPoint_Local, iVar, Sol[iVar+2*nVar]); } - if (fluid_structure && !dynamic) { + if (fluid_structure) { nodes->SetSolution_Pred(iPoint_Local, iVar, Sol[iVar]); nodes->SetSolution_Pred_Old(iPoint_Local, iVar, Sol[iVar]); + if (dynamic) nodes->SetSolution_Vel_Pred(iPoint_Local, iVar, Sol[iVar+nVar]); } if (fluid_structure && discrete_adjoint){ nodes->SetSolution_Old(iPoint_Local, iVar, Sol[iVar]); @@ -3331,12 +3369,16 @@ void CFEASolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig *c solver[MESH_0][FEA_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION_FEA_OLD); solver[MESH_0][FEA_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION_FEA_OLD); } - if (fluid_structure && !dynamic) { + if (fluid_structure) { solver[MESH_0][FEA_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION_PRED); solver[MESH_0][FEA_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION_PRED); solver[MESH_0][FEA_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION_PRED_OLD); solver[MESH_0][FEA_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION_PRED_OLD); + if (dynamic) { + solver[MESH_0][FEA_SOL]->InitiateComms(geometry[MESH_0], config, SOLUTION_VEL_PRED); + solver[MESH_0][FEA_SOL]->CompleteComms(geometry[MESH_0], config, SOLUTION_VEL_PRED); + } } /*--- Delete the class memory that is used to load the restart. ---*/ diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 2c2bd1ad4f50..291569a17950 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -513,7 +513,10 @@ void CMeshSolver::DeformMesh(CGeometry **geometry, CNumerics **numerics, CConfig UpdateDualGrid(geometry[MESH_0], config); /*--- The Grid Velocity is only computed if the problem is time domain ---*/ - if (time_domain) ComputeGridVelocity(geometry[MESH_0], config); + if (time_domain) { + if (config->GetFSI_Simulation()) ComputeGridVelocity_FromBoundary(geometry, numerics, config); + else ComputeGridVelocity(geometry[MESH_0], config); + } /*--- Update the multigrid structure. ---*/ UpdateMultiGrid(geometry, config); @@ -564,6 +567,43 @@ void CMeshSolver::UpdateDualGrid(CGeometry *geometry, CConfig *config){ } +void CMeshSolver::ComputeGridVelocity_FromBoundary(CGeometry **geometry, CNumerics **numerics, CConfig *config){ + + /*--- Compute the stiffness matrix, no point recording because we clear the residual. ---*/ + + const bool wasActive = AD::BeginPassive(); + + Compute_StiffMatrix(geometry[MESH_0], numerics, config); + + AD::EndPassive(wasActive); + + /*--- Clear residual (loses AD info), we do not want an incremental solution. ---*/ + SU2_OMP_PARALLEL { + LinSysRes.SetValZero(); + LinSysSol.SetValZero(); + } + + /*--- Impose boundary conditions including bouundary velocity ---*/ + SetBoundaryVelocities(geometry[MESH_0], numerics[FEA_TERM], config); + + /*--- Solve the linear system. ---*/ + Solve_System(geometry[MESH_0], config); + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + su2double val_vel = LinSysSol(iPoint, iDim); + + /*--- Non-dimensionalize velocity ---*/ + val_vel = val_vel/config->GetVelocity_Ref(); + + geometry[MESH_0]->nodes->SetGridVel(iPoint, iDim, val_vel); + } + } + + /*--- The velocity was computed for nPointDomain, now we communicate it. ---*/ + geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, GRID_VELOCITY); + geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, GRID_VELOCITY); +} + void CMeshSolver::ComputeGridVelocity(CGeometry *geometry, CConfig *config){ /*--- Compute the velocity of each node in the domain of the current rank @@ -707,6 +747,59 @@ void CMeshSolver::SetBoundaryDisplacements(CGeometry *geometry, CNumerics *numer } +void CMeshSolver::SetBoundaryVelocities(CGeometry *geometry, CNumerics *numerics, CConfig *config){ + + unsigned short iMarker; + + /*--- Impose zero displacements of all non-moving surfaces (also at nodes in multiple moving/non-moving boundaries). ---*/ + /*--- 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_Moving(iMarker) == NO) && + (config->GetMarker_All_KindBC(iMarker) != SYMMETRY_PLANE) && + (config->GetMarker_All_KindBC(iMarker) != SEND_RECEIVE) && + (config->GetMarker_All_KindBC(iMarker) != PERIODIC_BOUNDARY)) { + + BC_Clamped(geometry, numerics, config, iMarker); + } + } + + /*--- 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_Moving(iMarker) == NO) && + (config->GetMarker_All_KindBC(iMarker) == SYMMETRY_PLANE)) { + + BC_Clamped(geometry, numerics, config, iMarker); + } + } + + /*--- Impose velocity boundary conditions. ---*/ + for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { + if ((config->GetMarker_All_Deform_Mesh(iMarker) == YES) || + (config->GetMarker_All_Moving(iMarker) == YES)) { + + BC_Velocity(geometry, numerics, config, iMarker); + } + } + + /*--- Clamp far away nodes according to deform limit. ---*/ + if ((config->GetDeform_Stiffness_Type() == SOLID_WALL_DISTANCE) && + (config->GetDeform_Limit() < MaxDistance)) { + + const su2double limit = config->GetDeform_Limit() / MaxDistance; + + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { + if (nodes->GetWallDistance(iPoint) <= limit) continue; + + su2double zeros[MAXNVAR] = {0.0}; + nodes->SetSolution(iPoint, zeros); + LinSysSol.SetBlock(iPoint, zeros); + Jacobian.EnforceSolutionAtNode(iPoint, zeros, LinSysRes); + } + } +} + void CMeshSolver::SetDualTime_Mesh(void){ nodes->Set_Solution_time_n1(); @@ -754,6 +847,8 @@ void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * minus the coordinates of the reference mesh file ---*/ su2double displ = curr_coord - nodes->GetMesh_Coord(iPoint_Local, iDim); nodes->SetSolution(iPoint_Local, iDim, displ); + su2double vel = Restart_Data[index+iDim+6]; + if (time_domain && config->GetFSI_Simulation()) geometry[MESH_0]->nodes->SetGridVel(iPoint_Local, iDim, Restart_Data[index+iDim+6]); } /*--- Increment the overall counter for how many points have been loaded. ---*/ @@ -789,7 +884,7 @@ void CMeshSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig * UpdateDualGrid(geometry[MESH_0], config); /*--- For time-domain problems, we need to compute the grid velocities ---*/ - if (time_domain){ + if (time_domain && !config->GetFSI_Simulation()){ /*--- Update the old geometry (coordinates n and n-1) ---*/ Restart_OldGeometry(geometry[MESH_0], config); /*--- Once Displacement_n and Displacement_n1 are filled, @@ -837,7 +932,7 @@ void CMeshSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config) { /*--- Multizone problems require the number of the zone to be appended. ---*/ - if (nZone > 1) + if (config->GetMultizone_Problem()) filename = config->GetMultizone_FileName(filename, iZone, ""); /*--- Determine how many files need to be read. ---*/ @@ -855,11 +950,24 @@ void CMeshSolver::Restart_OldGeometry(CGeometry *geometry, CConfig *config) { if (Unst_RestartIter < 0) { - if (rank == MASTER_NODE) cout << "Requested mesh restart filename is negative. Setting known solution" << endl; + if (rank == MASTER_NODE) cout << "Requested mesh restart filename is negative. Setting zero displacement" << endl; + + unsigned long iPoint_Global; + + for (iPoint_Global = 0; iPoint_Global < geometry->GetGlobal_nPointDomain(); iPoint_Global++) { - /*--- Set loaded solution into correct previous time containers. ---*/ - if(iStep==1) nodes->Set_Solution_time_n(); - else nodes->Set_Solution_time_n1(); + auto iPoint_Local = geometry->GetGlobal_to_Local_Point(iPoint_Global); + + if (iPoint_Local >= 0) { + + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + if(iStep==1) + nodes->Set_Solution_time_n(iPoint_Local, iDim, 0.0); + else + nodes->Set_Solution_time_n1(iPoint_Local, iDim, 0.0); + } + } + } } else { string filename_n = config->GetUnsteady_FileName(filename, Unst_RestartIter, ""); diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 17393acf470c..7e25990344e3 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -1636,12 +1636,19 @@ void CSolver::GetCommCountAndType(const CConfig* config, COUNT_PER_POINT = nVar*3; MPI_TYPE = COMM_TYPE_DOUBLE; break; + case SOLUTION_VEL_PRED: + COUNT_PER_POINT = nVar; + MPI_TYPE = COMM_TYPE_DOUBLE; + break; case SOLUTION_PRED: COUNT_PER_POINT = nVar; MPI_TYPE = COMM_TYPE_DOUBLE; break; case SOLUTION_PRED_OLD: - COUNT_PER_POINT = nVar*3; + if (config->GetTime_Domain()) + COUNT_PER_POINT = nVar*4; + else + COUNT_PER_POINT = nVar*3; MPI_TYPE = COMM_TYPE_DOUBLE; break; case AUXVAR_GRADIENT: @@ -1798,6 +1805,10 @@ void CSolver::InitiateComms(CGeometry *geometry, bufDSend[buf_offset+nVar*2+iVar] = base_nodes->GetSolution_Accel_time_n(iPoint, iVar); } break; + case SOLUTION_VEL_PRED: + for (iVar = 0; iVar < nVar; iVar++) + bufDSend[buf_offset+iVar] = base_nodes->GetSolution_Vel_Pred(iPoint, iVar); + break; case SOLUTION_PRED: for (iVar = 0; iVar < nVar; iVar++) bufDSend[buf_offset+iVar] = base_nodes->GetSolution_Pred(iPoint, iVar); @@ -1807,6 +1818,8 @@ void CSolver::InitiateComms(CGeometry *geometry, bufDSend[buf_offset+iVar] = base_nodes->GetSolution_Old(iPoint, iVar); bufDSend[buf_offset+nVar+iVar] = base_nodes->GetSolution_Pred(iPoint, iVar); bufDSend[buf_offset+nVar*2+iVar] = base_nodes->GetSolution_Pred_Old(iPoint, iVar); + if (config->GetTime_Domain()) + bufDSend[buf_offset+nVar*3+iVar] = base_nodes->GetSolution_Vel_Pred(iPoint, iVar); } break; case MESH_DISPLACEMENTS: @@ -1976,6 +1989,10 @@ void CSolver::CompleteComms(CGeometry *geometry, base_nodes->SetSolution_Accel_time_n(iPoint, iVar, bufDRecv[buf_offset+nVar*2+iVar]); } break; + case SOLUTION_VEL_PRED: + for (iVar = 0; iVar < nVar; iVar++) + base_nodes->SetSolution_Vel_Pred(iPoint, iVar, bufDRecv[buf_offset+iVar]); + break; case SOLUTION_PRED: for (iVar = 0; iVar < nVar; iVar++) base_nodes->SetSolution_Pred(iPoint, iVar, bufDRecv[buf_offset+iVar]); @@ -1985,6 +2002,8 @@ void CSolver::CompleteComms(CGeometry *geometry, base_nodes->SetSolution_Old(iPoint, iVar, bufDRecv[buf_offset+iVar]); base_nodes->SetSolution_Pred(iPoint, iVar, bufDRecv[buf_offset+nVar+iVar]); base_nodes->SetSolution_Pred_Old(iPoint, iVar, bufDRecv[buf_offset+nVar*2+iVar]); + if (config->GetTime_Domain()) + base_nodes->SetSolution_Vel_Pred(iPoint, iVar, bufDRecv[buf_offset+nVar*3+iVar]); } break; case MESH_DISPLACEMENTS: diff --git a/SU2_CFD/src/variables/CDiscAdjFEABoundVariable.cpp b/SU2_CFD/src/variables/CDiscAdjFEABoundVariable.cpp index 4e929d4a0fdc..a442d0fae687 100644 --- a/SU2_CFD/src/variables/CDiscAdjFEABoundVariable.cpp +++ b/SU2_CFD/src/variables/CDiscAdjFEABoundVariable.cpp @@ -47,5 +47,6 @@ void CDiscAdjFEABoundVariable::AllocateBoundaryVariables(CConfig *config) { FlowTraction_Sens.resize(nBoundPt,nDim) = su2double(0.0); SourceTerm_DispAdjoint.resize(nBoundPt,nDim) = su2double(0.0); + SourceTerm_VelAdjoint.resize(nBoundPt,nDim) = su2double(0.0); } diff --git a/SU2_CFD/src/variables/CDiscAdjFEAVariable.cpp b/SU2_CFD/src/variables/CDiscAdjFEAVariable.cpp index ea6a9d708527..792c9919f22b 100644 --- a/SU2_CFD/src/variables/CDiscAdjFEAVariable.cpp +++ b/SU2_CFD/src/variables/CDiscAdjFEAVariable.cpp @@ -35,6 +35,7 @@ CDiscAdjFEAVariable::CDiscAdjFEAVariable(const su2double *disp, const su2double Solution_Direct.resize(nPoint,nVar); Sensitivity.resize(nPoint,nDim) = su2double(0.0); + Sensitivity_Old.resize(nPoint,nDim) = su2double(0.0); for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) for (unsigned long iVar = 0; iVar < nVar; iVar++) diff --git a/SU2_CFD/src/variables/CDiscAdjVariable.cpp b/SU2_CFD/src/variables/CDiscAdjVariable.cpp index 25eaf6433c66..bacaffa43257 100644 --- a/SU2_CFD/src/variables/CDiscAdjVariable.cpp +++ b/SU2_CFD/src/variables/CDiscAdjVariable.cpp @@ -47,6 +47,7 @@ CDiscAdjVariable::CDiscAdjVariable(const su2double* sol, unsigned long npoint, u Solution_Direct.resize(nPoint,nVar); Sensitivity.resize(nPoint,nDim) = su2double(0.0); + Sensitivity_Old.resize(nPoint,nDim) = su2double(0.0); for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) for (unsigned long iVar = 0; iVar < nVar; ++iVar) diff --git a/SU2_CFD/src/variables/CFEAVariable.cpp b/SU2_CFD/src/variables/CFEAVariable.cpp index 376db841bc69..b3f0971da340 100644 --- a/SU2_CFD/src/variables/CFEAVariable.cpp +++ b/SU2_CFD/src/variables/CFEAVariable.cpp @@ -64,11 +64,23 @@ CFEAVariable::CFEAVariable(const su2double *val_fea, unsigned long npoint, unsig } Solution_Vel_time_n = Solution_Vel; Solution_Accel_time_n = Solution_Accel; + + if(config->GetMultizone_Problem() && config->GetAD_Mode()) { + AD_Vel_InputIndex.resize(nPoint,nVar) = -1; + AD_Vel_OutputIndex.resize(nPoint,nVar) = -1; + AD_Vel_Time_n_InputIndex.resize(nPoint,nVar) = -1; + AD_Vel_Time_n_OutputIndex.resize(nPoint,nVar) = -1; + AD_Accel_InputIndex.resize(nPoint,nVar) = -1; + AD_Accel_OutputIndex.resize(nPoint,nVar) = -1; + AD_Accel_Time_n_InputIndex.resize(nPoint,nVar) = -1; + AD_Accel_Time_n_OutputIndex.resize(nPoint,nVar) = -1; + } } if (fsi_analysis) { Solution_Pred = Solution; Solution_Pred_Old = Solution; + if (dynamic_analysis) Solution_Vel_Pred = Solution_Vel; } /*--- If we are going to use incremental analysis, we need a way to store the old solution ---*/ @@ -98,46 +110,107 @@ void CFEAVariable::SetSolution_Vel_time_n() { Solution_Vel_time_n = Solution_Vel void CFEAVariable::SetSolution_Accel_time_n() { Solution_Accel_time_n = Solution_Accel; } -void CFEAVariable::Register_femSolution_time_n() { - for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) - for (unsigned long iVar = 0; iVar < nVar; iVar++) - AD::RegisterInput(Solution_time_n(iPoint,iVar)); +void CFEAVariable::Register_femSolution_time_n(bool input, bool push_index) { + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + for(unsigned long iVar=0; iVarGetTime_Domain()) Boundary_Velocity.resize(nBoundPt,nDim) = su2double(0.0); } void CMeshBoundVariable::Register_BoundDisp(bool input) { diff --git a/SU2_CFD/src/variables/CVariable.cpp b/SU2_CFD/src/variables/CVariable.cpp index 917d4f3fed28..5a09c2c15609 100644 --- a/SU2_CFD/src/variables/CVariable.cpp +++ b/SU2_CFD/src/variables/CVariable.cpp @@ -74,6 +74,12 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva if(config->GetMultizone_Problem() && config->GetAD_Mode()) { AD_InputIndex.resize(nPoint,nVar) = -1; AD_OutputIndex.resize(nPoint,nVar) = -1; + if (config->GetTime_Domain()) { + AD_Time_n_InputIndex.resize(nPoint,nVar) = -1; + AD_Time_n_OutputIndex.resize(nPoint,nVar) = -1; + AD_Time_n1_InputIndex.resize(nPoint,nVar) = -1; + AD_Time_n1_OutputIndex.resize(nPoint,nVar) = -1; + } } if (config->GetMultizone_Problem()) @@ -135,14 +141,30 @@ void CVariable::RegisterSolution(bool input, bool push_index) { } } -void CVariable::RegisterSolution_time_n() { - for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) - for(unsigned long iVar=0; iVar