diff --git a/AUTHORS.md b/AUTHORS.md index f6580fadb89..c530c7dc9f0 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -82,6 +82,7 @@ Giacomo Baldan Giulio Gori Guillaume Bâty Harichand M V +Harsh Patel HL Kline HL Zhi IndianaStokes diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 7405a527eac..89f56de6d68 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1095,8 +1095,9 @@ class CConfig { long ParMETIS_pointWgt; /*!< \brief Load balancing weight given to points. */ long ParMETIS_edgeWgt; /*!< \brief Load balancing weight given to edges. */ unsigned short DirectDiff; /*!< \brief Direct Differentation mode. */ - bool DiscreteAdjoint, /*!< \brief AD-based discrete adjoint mode. */ - DiscreteAdjointDebug; /*!< \brief Discrete adjoint debug mode using tags. */ + bool DiscreteAdjoint; /*!< \brief AD-based discrete adjoint mode. */ + ENUM_DISC_ADJ_TYPE Kind_DiscreteAdjoint; /*!< \brief AD-based discrete adjoint formulation. */ + bool DiscreteAdjointDebug; /*!< \brief Discrete adjoint debug mode using tags. */ su2double Const_DES; /*!< \brief Detached Eddy Simulation Constant. */ WINDOW_FUNCTION Kind_WindowFct; /*!< \brief Type of window (weight) function for objective functional. */ unsigned short Kind_HybridRANSLES; /*!< \brief Kind of Hybrid RANS/LES. */ @@ -9188,6 +9189,12 @@ class CConfig { */ bool GetDiscrete_Adjoint(void) const { return DiscreteAdjoint; } + /*! + * \brief Get the kind of discrete adjoint solver formulation. + * \return the discrete adjoint indicator. + */ + ENUM_DISC_ADJ_TYPE GetKind_DiscreteAdjoint(void) const { return Kind_DiscreteAdjoint; } + /*! * \brief Get the indicator whether a debug run for the discrete adjoint solver will be started. * \return the discrete adjoint debug indicator. diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index 563674def66..e45418f8129 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -1400,6 +1400,16 @@ class CGeometry { */ inline virtual void SetSensitivity(unsigned long iPoint, unsigned short iDim, su2double val) {} + /*! + * \brief Get a reference to the sensitivity matrix. + * \return Reference to the sensitivity matrix. + */ + inline virtual const su2activematrix& GetSensitivityMatrix() const { + SU2_MPI::Error("GetSensitivityMatrix not implemented for this geometry type.", CURRENT_FUNCTION); + static su2activematrix sens_matrix; + return sens_matrix; + } + /*! * \brief Get the average normal at a specific span for a given marker in the turbomachinery reference of frame. * \param[in] val_marker - marker value. diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp index be0863584d4..ae75737176a 100644 --- a/Common/include/geometry/CPhysicalGeometry.hpp +++ b/Common/include/geometry/CPhysicalGeometry.hpp @@ -779,6 +779,12 @@ class CPhysicalGeometry final : public CGeometry { Sensitivity(iPoint, iDim) = val; } + /*! + * \brief Get a reference to the sensitivity matrix. + * \return Reference to the sensitivity matrix. + */ + inline const su2activematrix& GetSensitivityMatrix() const override { return Sensitivity; } + /*! * \brief Check the mesh for periodicity and deactivate multigrid if periodicity is found. * \param[in] config - Definition of the particular problem. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 1361a12072d..ad8f8075fd0 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -2798,6 +2798,25 @@ static const MapType Streamwise_Periodic_ MakePair("MASSFLOW", ENUM_STREAMWISE_PERIODIC::MASSFLOW) }; +/*! + * \brief Types of discrete adjoint solver formulations. + */ +enum class ENUM_DISC_ADJ_TYPE { + FIXED_POINT, /*!< \brief Fixed-point discrete-adjoint formulation. */ + RESIDUALS /*!< \brief Residual-based discrete-adjoint formulation. */ +}; +static const MapType DiscreteAdjoint_Map = { + MakePair("FIXED_POINT", ENUM_DISC_ADJ_TYPE::FIXED_POINT) + MakePair("RESIDUALS", ENUM_DISC_ADJ_TYPE::RESIDUALS) +}; + +enum class ENUM_VARIABLE { + RESIDUALS, + OBJECTIVE, + TRACTIONS, + COORDINATES +}; + /*! * \brief Container to hold Variables for streamwise Periodic flow as they are often used together in places. */ diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 79fe6d6aab1..3007286db82 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1140,6 +1140,8 @@ void CConfig::SetConfig_Options() { #endif /*!\brief MATH_PROBLEM \n DESCRIPTION: Mathematical problem \n Options: DIRECT, ADJOINT \ingroup Config*/ addMathProblemOption("MATH_PROBLEM", ContinuousAdjoint, false, DiscreteAdjoint, discAdjDefault, Restart_Flow, discAdjDefault); + /*!\brief DISC_ADJ_TYPE \n DESCRIPTION: Discrete adjoint problem formulation \n Options: FIXED_POINT, RESIDUALS \ingroup Config*/ + addEnumOption("KIND_DISC_ADJ", Kind_DiscreteAdjoint, DiscreteAdjoint_Map, ENUM_DISC_ADJ_TYPE::FIXED_POINT); /*!\brief KIND_TURB_MODEL \n DESCRIPTION: Specify turbulence model \n Options: see \link Turb_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/ addEnumOption("KIND_TURB_MODEL", Kind_Turb_Model, Turb_Model_Map, TURB_MODEL::NONE); /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SST turbulence model options/corrections. \n Options: see \link SST_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/ diff --git a/SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp b/SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp index 7d7e8aa8d3e..41af6795f7f 100644 --- a/SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp +++ b/SU2_CFD/include/drivers/CDiscAdjSinglezoneDriver.hpp @@ -2,7 +2,7 @@ * \file CDiscAdjSinglezoneDriver.hpp * \brief Headers of the main subroutines for driving single or multi-zone problems. * The subroutines and functions are in the driver_structure.cpp file. - * \author T. Economon, H. Kline, R. Sanchez + * \author T. Economon, H. Kline, R. Sanchez, H. Patel, A. Gastaldi * \version 8.4.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -28,6 +28,19 @@ #pragma once #include "CSinglezoneDriver.hpp" +#include "../../../Common/include/toolboxes/CQuasiNewtonInvLeastSquares.hpp" +#include "../../../Common/include/linear_algebra/CPreconditioner.hpp" +#include "../../../Common/include/linear_algebra/CMatrixVectorProduct.hpp" +#include "../../../Common/include/linear_algebra/CSysSolve.hpp" + +#ifdef CODI_FORWARD_TYPE + using Scalar = su2double; +#else + using Scalar = su2mixedfloat; +#endif + +class LinOperator; +class LinPreconditioner; /*! * \class CDiscAdjSinglezoneDriver @@ -55,6 +68,24 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver { COutput *direct_output; CNumerics ***numerics; /*!< \brief Container vector with all the numerics. */ + /*!< \brief Fixed-Point corrector that can be applied to inner iterations of the residual-based formulation. */ + CQuasiNewtonInvLeastSquares Corrector; + + /*!< \brief Members to use GMRES to drive inner iterations (alternative to quasi-Newton) of the residual-based formulation. */ + static constexpr unsigned long KrylovMinIters = 3; + const Scalar KrylovSysTol = 0.00001; + const Scalar KrylovPreTol = 0.1; + bool KrylovSet = false; + + CSysMatrix CopiedJacobian; + CSysSolve AdjSolver; + CSysVector AdjRHS; + CSysVector AdjSol; + CPreconditioner* PrimalPreconditioner = nullptr; + CSysMatrixVectorProduct* PrimalJacobian = nullptr; + LinOperator* AdjOperator = nullptr; + LinPreconditioner* AdjPreconditioner = nullptr; + /*! * \brief Record one iteration of a flow iteration in within multiple zones. * \param[in] kind_recording - Type of recording (full list in ENUM_RECORDING, option_structure.hpp) @@ -65,27 +96,58 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver { * \brief Run one iteration of the solver. * \param[in] kind_recording - Type of recording (full list in ENUM_RECORDING, option_structure.hpp) */ - void DirectRun(RECORDING kind_recording); + void DirectRunFixedPoint(RECORDING kind_recording); + + /*! + * \brief Run one iteration of the solver. + * \param[in] kind_recording - Type of recording (full list in ENUM_RECORDING, option_structure.hpp) + */ + void DirectRunResidual(RECORDING kind_recording); + + /*! + * \brief Run a single iteration of the main fixed-point discrete adjoint solver with a single zone. + */ + void RunFixedPoint(); + + /*! + * \brief Run the computation of the main residual-based discrete adjoint sensitivities with a single zone. + */ + void RunResidual(); + + /*! + * \brief Run a single iteration of the secondary fixed-point discrete adjoint solver with a single zone. + */ + void SecondaryRunFixedPoint(); + + /*! + * \brief Run the computation of the secondary residual-based discrete adjoint sensitivities with a single zone. + */ + void SecondaryRunResidual(); /*! - * \brief Set the objective function. + * \brief Update the fixed-point discrete adjoint solver with a single zone. */ - void SetObjFunction(void); + void UpdateAdjointsFixedPoint(); + + /*! + * \brief Update the residual-based discrete adjoint solver with a single zone. + */ + void UpdateAdjointsResidual(); /*! * \brief Initialize the adjoint value of the objective function. */ - void SetAdjObjFunction(void); + void SetAdjointObjective(); /*! * \brief Record the main computational path. */ - void MainRecording(void); + void MainRecording(); /*! * \brief Record the secondary computational path. */ - void SecondaryRecording(void); + void SecondaryRecording(); /*! * \brief gets Convergence on physical time scale, (deactivated in adjoint case) @@ -99,7 +161,6 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver { * \brief Constructor of the class. * \param[in] confFile - Configuration file name. * \param[in] val_nZone - Total number of zones. - * \param[in] val_nDim - Total number of dimensions. * \param[in] MPICommunicator - MPI communicator for SU2. */ CDiscAdjSinglezoneDriver(char* confFile, @@ -120,10 +181,149 @@ class CDiscAdjSinglezoneDriver : public CSinglezoneDriver { /*! * \brief Run a single iteration of the discrete adjoint solver with a single zone. */ - void Run(void) override; + void Run() override; + + /*! + * \brief Update the discrete adjoint solution. + */ + void UpdateAdjoints(); + + /*! + * \brief Update the primal time iteration. + */ + void UpdateTimeIter(); + + /*! + * \brief Update the primal far-field variables. + */ + void UpdateFarfield(); + + /*! + * \brief Update the primal geometry (does not include mesh deformation). + */ + void UpdateGeometry(); + + /*! + * \brief Deform the primal mesh. + */ + void DeformGeometry(); + + /*! + * \brief Update the primal states. + */ + void UpdateStates(); + + /*! + * \brief Update the primal residuals. + */ + void UpdateResiduals(); + + /*! + * \brief Set wall-normal constraint R_n = (rho*v).n so that (dR/dq)^T is non-singular. + */ + void SetWallNormalConstraint(); + + /*! + * \brief Update the primal tractions. + */ + void UpdateTractions(); + + /*! + * \brief Update the primal objective. + */ + void UpdateObjective(); + + /*! + * \brief Update the primal jacobian. + */ + void UpdateJacobians(); /*! * \brief Postprocess the adjoint iteration for ZONE_0. */ - void Postprocess(void) override; + void Postprocess() override; + + /*! + * \brief Get the partial objective sensitivities from all solvers. + * \param[out] values - Values object with interface (iPoint, iVar). + */ + template + void GetAllObjectiveStatesSensitivities(Container& values) const { + const auto nPoint = geometry_container[ZONE_0][INST_0][MESH_0]->GetnPoint(); + + /*--- Get all the partial sensitivities (dObjective/dStates) ---*/ + unsigned short offset = 0; + + for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) { + auto solver = solver_container[ZONE_0][INST_0][MESH_0][iSol]; + + if (!solver || !solver->GetAdjoint()) continue; + + const auto& mat = solver->GetPartialMatrix_dObjective_dStates(); + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { + for (auto iVar = 0ul; iVar < solver->GetnVar(); ++iVar) { + values(iPoint, offset + iVar) = -SU2_TYPE::GetValue(mat(iPoint, iVar)); + } + } + + offset += solver->GetnVar(); + } + } + + /*! + * \brief Get the partial jacobian-adjoint products from all solvers. + * \param[out] values - Values object with interface (iPoint, iVar). + */ + template + void GetAllResidualsStatesSensitivities(Container& values) const { + const auto nPoint = geometry_container[ZONE_0][INST_0][MESH_0]->GetnPoint(); + + /*--- Get all the partial jacobian-adjoint products (dResiduals/dStates) ---*/ + unsigned short offset = 0; + + for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) { + auto solver = solver_container[ZONE_0][INST_0][MESH_0][iSol]; + + if (!solver || !solver->GetAdjoint()) continue; + + const auto& mat = solver->GetPartialMatrix_dResiduals_dStates(); + for (auto iPoint = 0ul; iPoint < nPoint; ++iPoint) { + for (auto iVar = 0ul; iVar < solver->GetnVar(); ++iVar) { + values(iPoint, offset + iVar) = SU2_TYPE::GetValue(mat(iPoint, iVar)); + } + } + + offset += solver->GetnVar(); + } + } + + /*! + * \brief Adjoint problem Jacobian-vector product. + */ + void ApplyOperator(const CSysVector& u, CSysVector& v); + + /*! + * \brief Adjoint problem preconditioner (based on the transpose approximate Jacobian of the primal problem). + */ + void ApplyPreconditioner(const CSysVector& u, CSysVector& v); +}; + +class LinOperator : public CMatrixVectorProduct { + public: + CDiscAdjSinglezoneDriver* const driver; + LinOperator(CDiscAdjSinglezoneDriver* d) : driver(d) { } + + inline void operator()(const CSysVector & u, CSysVector & v) const override { + driver->ApplyOperator(u, v); + } +}; + +class LinPreconditioner : public CPreconditioner { + public: + CDiscAdjSinglezoneDriver* const driver; + LinPreconditioner(CDiscAdjSinglezoneDriver* d) : driver(d) { } + + inline void operator()(const CSysVector & u, CSysVector & v) const override { + driver->ApplyPreconditioner(u, v); + } }; diff --git a/SU2_CFD/include/drivers/CDriver.hpp b/SU2_CFD/include/drivers/CDriver.hpp index b01f560fa28..015339b55de 100644 --- a/SU2_CFD/include/drivers/CDriver.hpp +++ b/SU2_CFD/include/drivers/CDriver.hpp @@ -469,6 +469,16 @@ class CDriver : public CDriverBase { */ void BoundaryConditionsUpdate(); + /*! + * \brief Update the geometry (i.e. dual grid). + */ + void UpdateGeometry(); + + /*! + * \brief Update the primal far-field variables. + */ + void UpdateFarfield(); + /*! * \brief Get the number of time iterations. * \return Number of time iterations. @@ -509,17 +519,181 @@ class CDriver : public CDriverBase { */ void SetInletAngle(unsigned short iMarker, passivedouble alpha); + /*! + * \brief Preprocess the inlets via file input for all solvers. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + + /*! + * \brief Get the free-stream Reynolds number. + * \return Free-stream Reynolds number. + */ + passivedouble GetReynoldsNumber() const; + + /*! + * \brief Get the free-stream Mach number. + * \return Free-stream Mach number. + */ + passivedouble GetMachNumber() const; + + /*! + * \brief Get the free-stream angle of attack (in degrees). + * \return Free-stream angle of attack. + */ + passivedouble GetAngleOfAttack() const; + + /*! + * \brief Get the free-stream angle of side-slip (in degrees). + * \return Free-stream angle of side-slip. + */ + passivedouble GetAngleOfSideslip() const; + + /*! + * \brief Set the free-stream Reynolds number. + * \param[in] value - User-defined Reynolds number. + */ + void SetReynoldsNumber(passivedouble value); + + /*! + * \brief Set the free-stream Mach number. + * \param[in] value - User-defined Mach number. + */ + void SetMachNumber(passivedouble value); + /*! * \brief Set the angle of attack of the farfield. * \param[in] alpha - Angle (degree). */ - void SetFarFieldAoA(passivedouble alpha); + void SetAngleOfAttack(passivedouble alpha); /*! * \brief Set the angle of sideslip of the farfield. * \param[in] beta - Angle (degree). */ - void SetFarFieldAoS(passivedouble beta); + void SetAngleOfSideslip(passivedouble beta); + + /*! + * \brief Get the number of conservative state variables. + * \return Number of conservative state variables. + */ + unsigned long GetNumberStateVariables() const; + + /*! + * \brief Get the number of primitive state variables. + * \return Number of primitive state variables. + */ + unsigned long GetNumberPrimitiveVariables() const; + + /*! + * \brief Get the adjoint flow forces at a marker vertex. + * \param[in] iMarker - Marker index. + * \return CPyWrapperMatrixView of shape (nVertex, nDim). + */ + CPyWrapperMatrixView MarkerAdjointForces(unsigned short iMarker) const; + + /*! + * \brief Get a read-only view of dCoordinates/dCoordinates^T * psi for all mesh nodes. + * \return CPyWrapperMatrixView of shape (nPoint, nDim). + */ + CPyWrapperMatrixView CoordinatesCoordinatesSensitivities() const; + + /*! + * \brief Get a read-only view of dCoordinates/dDisplacements^T * psi for a marker. + * \param[in] iMarker - Marker index. + * \return CPyWrapperMatrixView of shape (nVertex, nDim). + */ + CPyWrapperMatrixView MarkerCoordinatesDisplacementsSensitivities(unsigned short iMarker) const; + + /*! + * \brief Get sensitivity of objective function with respect to farfield design variables as a partial derivative. + * \return Partial derivative of aerodynamic function with respect to farfield design variable. + */ + vector GetObjectiveFarfieldVariablesSensitivities() const; + + /*! + * \brief Get sensitivity of flow residuals with respect to farfield design variables as a matrix-vector product with + * the adjoint variable. + * \return Partial derivative of aerodynamic residuals with respect to farfield design variable. + */ + vector GetResidualsFarfieldVariablesSensitivities() const; + + /*! + * \brief Get a read-only view of dObjective/dStates for all mesh nodes. + * \return CPyWrapperMatrixView of shape (nPoint, nVar). + */ + CPyWrapperMatrixView ObjectiveStatesSensitivities() const; + + /*! + * \brief Get a read-only view of dResiduals/dStates^T * psi for all mesh nodes. + * \return CPyWrapperMatrixView of shape (nPoint, nVar). + */ + CPyWrapperMatrixView ResidualsStatesSensitivities() const; + + /*! + * \brief Get a read-only view of dTractions/dStates^T * psi for all mesh nodes. + * \return CPyWrapperMatrixView of shape (nPoint, nVar). + */ + CPyWrapperMatrixView ForcesStatesSensitivities() const; + + /*! + * \brief Get a read-only view of dObjective/dCoordinates for all mesh nodes. + * \return CPyWrapperMatrixView of shape (nPoint, nDim). + */ + CPyWrapperMatrixView ObjectiveCoordinatesSensitivities() const; + + /*! + * \brief Get a read-only view of dResiduals/dCoordinates^T * psi for all mesh nodes. + * \return CPyWrapperMatrixView of shape (nPoint, nDim). + */ + CPyWrapperMatrixView ResidualsCoordinatesSensitivities() const; + + /*! + * \brief Get a read-only view of dTractions/dCoordinates^T * psi for all mesh nodes. + * \return CPyWrapperMatrixView of shape (nPoint, nDim). + */ + CPyWrapperMatrixView ForcesCoordinatesSensitivities() const; + + /*! + * \brief Get a read-only view of dObjective/dDisplacements for a marker. + * \param[in] iMarker - Marker index. + * \return CPyWrapperMatrixView of shape (nVertex, nDim). + */ + CPyWrapperMatrixView MarkerObjectiveDisplacementsSensitivities(unsigned short iMarker) const; + + /*! + * \brief Get a read-only view of dResiduals/dDisplacements^T * psi for a marker. + * \param[in] iMarker - Marker index. + * \return CPyWrapperMatrixView of shape (nVertex, nDim). + */ + CPyWrapperMatrixView MarkerResidualsDisplacementsSensitivities(unsigned short iMarker) const; + + /*! + * \brief Get a read-only view of dTractions/dDisplacements^T * psi for a marker. + * \param[in] iMarker - Marker index. + * \return CPyWrapperMatrixView of shape (nVertex, nDim). + */ + CPyWrapperMatrixView MarkerForcesDisplacementsSensitivities(unsigned short iMarker) const; + + /*! + * \brief Get sensitivities of the flow forces for the structural solver. + * \param[in] iMarker - Marker index. + * \return Sensitivity of flow forces (nVertex, nDim). + */ + vector GetMarkerForceSensitivities(unsigned short iMarker) const; + + /*! + * \brief Get a read/write view of the adjoint source term for all mesh nodes. + * \return CPyWrapperMatrixView of shape (nPoint, nVar). + */ + CPyWrapperMatrixView AdjointSourceTerm(); + + /*! + * \brief Get all the flow load boundary marker tags. + * \return List of flow load boundary markers tags. + */ + vector GetFluidLoadMarkerTags() const; /*! * \brief Set the dynamic mesh translation rates. diff --git a/SU2_CFD/include/integration/CIntegration.hpp b/SU2_CFD/include/integration/CIntegration.hpp index f6f35910e67..80a53895149 100644 --- a/SU2_CFD/include/integration/CIntegration.hpp +++ b/SU2_CFD/include/integration/CIntegration.hpp @@ -137,4 +137,16 @@ class CIntegration { CNumerics ******numerics_container, CConfig **config, unsigned short RunTime_EqSystem, unsigned short iZone, unsigned short iInst) { }; + /*! + * \brief A virtual member. + * \param[in] geometry_container - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] numerics_container - Description of the numerical method (the way in which the equations are solved). + * \param[in] config_container - Definition of the particular problem. + * \param[in] EqSystem - System of equations which is going to be solved. + * \param[in] iZone - Zone index. + * \param[in] iInst - Inst index. + */ + virtual void ComputeResiduals(CGeometry ****geometry_container, CSolver *****solvers_container, CNumerics ******numerics_container, + CConfig **config_container, unsigned short EqSystem, unsigned short iZone, unsigned short iInst); }; diff --git a/SU2_CFD/include/solvers/CDiscAdjMeshSolver.hpp b/SU2_CFD/include/solvers/CDiscAdjMeshSolver.hpp index 14d0ded3963..2d577d03b76 100644 --- a/SU2_CFD/include/solvers/CDiscAdjMeshSolver.hpp +++ b/SU2_CFD/include/solvers/CDiscAdjMeshSolver.hpp @@ -75,6 +75,13 @@ class CDiscAdjMeshSolver final : public CSolver { */ ~CDiscAdjMeshSolver() override; + /* + * \brief Register the mesh variables as outputs. + * \param[in] geometry_container - The geometry container holding all grid levels. + * \param[in] config_container - The particular config. + */ + void RegisterOutput(CGeometry *geometry, CConfig *config) override; + /*! * \brief Performs the preprocessing of the AD-based mesh adjoint solver. * Registers all necessary variables on the tape. Called while tape is active. @@ -83,6 +90,13 @@ class CDiscAdjMeshSolver final : public CSolver { */ void RegisterSolution(CGeometry *geometry, CConfig *config) override; + /*! + * \brief Seed the volume coordinate adjoints of the AD-based mesh adjoint solver. + * \param[in] geometry_container - The geometry container holding all grid levels. + * \param[in] config_container - The particular config. + */ + void SetAdjoint_Output(CGeometry *geometry, CConfig *config) override; + /*! * \brief Sets the adjoint values of the input variables of the flow (+turb.) iteration * after tape has been evaluated. diff --git a/SU2_CFD/include/solvers/CDiscAdjResidualSolver.hpp b/SU2_CFD/include/solvers/CDiscAdjResidualSolver.hpp new file mode 100755 index 00000000000..5ded89f80c4 --- /dev/null +++ b/SU2_CFD/include/solvers/CDiscAdjResidualSolver.hpp @@ -0,0 +1,308 @@ +/*! + * \file CDiscAdjResidualSolver.hpp + * \brief Headers of the residual-based adjoint solver class. + * \author H. Patel, A. Gastaldi + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include "../variables/CDiscAdjVariable.hpp" +#include "CSolver.hpp" + +class CDiscAdjResidualSolver final : public CSolver { + private: + static constexpr size_t MAXNDIM = 3; /*!< \brief Max number of space dimensions, used in some static arrays. */ + static constexpr size_t MAXNVAR = 32; /*!< \brief Max number of variables, for static arrays. */ + + static constexpr size_t OMP_MAX_SIZE = 1024; /*!< \brief Max chunk size for light point loops. */ + + unsigned long omp_chunk_size; /*!< \brief Chunk size used in light point loops. */ + + unsigned short KindDirect_Solver; + CSolver* direct_solver; + + const unsigned short nTrim = 2; + su2matrix Partial_Prod_dCoordinates_dCoordinates; /*!< \brief Partial sensitivity of the volume coordinates w.r.t. the boundary displacements (matrix-vector product with adjoint vector). */ + vector> Partial_Prod_dCoordinates_dDisplacements; /*!< \brief Partial sensitivity of the volume coordinates w.r.t. the boundary displacements (matrix-vector product with adjoint vector). */ + su2vector Partial_Sens_dObjective_dVariables; /*!< \brief Partial sensitivity of the objective w.r.t. the input variables. */ + su2vector Partial_Prod_dResiduals_dVariables; /*!< \brief Partial sensitivity of the residuals w.r.t. the input variables (matrix-vector product with adjoint vector). */ + su2matrix Partial_Sens_dObjective_dStates; /*!< \brief Partial sensitivity of the objective w.r.t. the solution. */ + su2matrix Partial_Prod_dResiduals_dStates; /*!< \brief Partial sensitivity of the residuals w.r.t. the solution (matrix-vector product with adjoint vector). */ + su2matrix Partial_Prod_dTractions_dStates; /*!< \brief Partial sensitivity of the surface tractions w.r.t. the solution (matrix-vector product with adjoint vector). */ + su2matrix Partial_Sens_dObjective_dCoordinates; /*!< \brief Partial sensitivity of the objective w.r.t. the boundary displacements . */ + su2matrix Partial_Prod_dResiduals_dCoordinates; /*!< \brief Partial sensitivity of the residuals w.r.t. the boundary displacements (matrix-vector product with adjoint vector). */ + su2matrix Partial_Prod_dTractions_dCoordinates; /*!< \brief Partial sensitivity of the tractions w.r.t. the boundary displacements (matrix-vector product with traction adjoints). */ + vector> Partial_Sens_dObjective_dDisplacements; /*!< \brief Partial sensitivity of the objective w.r.t. the boundary displacements . */ + vector> Partial_Prod_dResiduals_dDisplacements; /*!< \brief Partial sensitivity of the residuals w.r.t. the boundary displacements (matrix-vector product with adjoint vector). */ + vector> Partial_Prod_dTractions_dDisplacements; /*!< \brief Partial sensitivity of the tractions w.r.t. the boundary displacements (matrix-vector product with traction adjoints). */ + su2matrix AD_ResidualIndex; /*!< \brief Indices of Residual variables in the adjoint vector. */ + + vector> CSensitivity; /*!< \brief Shape sensitivity coefficient for each boundary and vertex. */ + vector Sens_Geo; /*!< \brief Total shape sensitivity for each monitored boundary. */ + su2double Total_Sens_Mach; /*!< \brief Total mach sensitivity coefficient for all the boundaries. */ + su2double Total_Sens_AoA; /*!< \brief Total angle of attack sensitivity coefficient for all the boundaries. */ + su2double Total_Sens_Geo; /*!< \brief Total shape sensitivity coefficient for all the boundaries. */ + su2double Total_Sens_Press; /*!< \brief Total farfield sensitivity to pressure. */ + su2double Total_Sens_Temp; /*!< \brief Total farfield sensitivity to temperature. */ + su2double Total_Sens_BPress; /*!< \brief Total sensitivity to outlet pressure. */ + su2double Total_Sens_Density; /*!< \brief Total sensitivity to initial density (incompressible). */ + su2double Total_Sens_ModVel; /*!< \brief Total sensitivity to inlet velocity (incompressible). */ + su2double Mach, Alpha, Beta, Pressure, Temperature, BPressure, ModVel; + su2double TemperatureRad, Total_Sens_Temp_Rad; + + CDiscAdjVariable* nodes = nullptr; /*!< \brief Highest level in variable hierarchy this solver can safely use. */ + + /*! + * \brief Return nodes to allow CSolver::base_nodes to be set. + */ + inline CVariable* GetBaseClassPointerToNodes() override { return nodes; } + + public: + /*! + * \brief Constructor of the class. + */ + CDiscAdjResidualSolver() = default; + + /*! + * \overload + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] solver - Initialize the discrete adjoint solver with the corresponding direct solver. + * \param[in] Kind_Solver - The kind of direct solver. + * \param[in] iMesh - Index of the mesh in multi-grid computations. + */ + CDiscAdjResidualSolver(CGeometry* geometry, CConfig* config, CSolver* solver, unsigned short Kind_Solver, + unsigned short iMesh); + + /*! + * \brief Destructor of the class. + */ + ~CDiscAdjResidualSolver() override; + + /*! + * \brief Performs the preprocessing of the adjoint AD-based solver. Registers all necessary variables on the tape. Called while tape is active. + * \param[in] geometry_container - The geometry container holding all grid levels. + * \param[in] config_container - The particular config. + */ + void RegisterSolution(CGeometry* geometry, CConfig* config) override; + + /*! + * \brief Performs the preprocessing of the adjoint AD-based solver. Registers all necessary variables that are + * output variables on the tape. Called while tape is active. + * \param[in] geometry_container - The geometry container holding all grid levels. + * \param[in] config_container - The particular config. + */ + void RegisterOutput(CGeometry* geometry, CConfig* config) override; + + /*! + * \brief Sets the adjoint values of the output of the flow (+turbulence) iteration before evaluation of the tape. + * \param[in] geometry - The geometrical definition of the problem. + * \param[in] config - The particular config. + */ + void SetAdjoint_Output(CGeometry* geometry, CConfig* config) override; + + /*! + * \brief A virtual member. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] output - Kind of output variables. + */ + void ExtractAdjoint_Solution(CGeometry* geometry, CConfig* config, ENUM_VARIABLE variable); + + /*! + * \brief Set the surface sensitivity. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetSurface_Sensitivity(CGeometry* geometry, CConfig* config) override; + + /*! + * \brief Extract and set the geometrical sensitivity. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + */ + void SetSensitivity(CGeometry* geometry, CConfig* config, CSolver*) override; + + /*! + * \brief Provide the total shape sensitivity coefficient. + * \return Value of the geometrical sensitivity coefficient (inviscid + viscous contribution). + */ + inline su2double GetTotal_Sens_Geo() const override { return Total_Sens_Geo; } + + /*! + * \brief Set the total Mach number sensitivity coefficient. + * \return Value of the Mach sensitivity coefficient (inviscid + viscous contribution). + */ + inline su2double GetTotal_Sens_Mach() const override { return Total_Sens_Mach; } + + /*! + * \brief Set the total angle of attack sensitivity coefficient. + * \return Value of the angle of attack sensitivity coefficient (inviscid + viscous contribution). + */ + inline su2double GetTotal_Sens_AoA() const override { return Total_Sens_AoA; } + + /*! + * \brief Set the total farfield pressure sensitivity coefficient. + * \return Value of the farfield pressure sensitivity coefficient (inviscid + viscous contribution). + */ + inline su2double GetTotal_Sens_Press() const override { return Total_Sens_Press; } + + /*! + * \brief Set the total farfield temperature sensitivity coefficient. + * \return Value of the farfield temperature sensitivity coefficient (inviscid + viscous contribution). + */ + inline su2double GetTotal_Sens_Temp() const override { return Total_Sens_Temp; } + + /*! + * \brief Get the total Back pressure number sensitivity coefficient. + * \return Value of the Back sensitivity coefficient (inviscid + viscous contribution). + */ + inline su2double GetTotal_Sens_BPress() const override { return Total_Sens_BPress; } + + /*! + * \brief Get the total density sensitivity coefficient. + * \return Value of the density sensitivity. + */ + inline su2double GetTotal_Sens_Density() const override { return Total_Sens_Density; } + + /*! + * \brief Get the total velocity magnitude sensitivity coefficient. + * \return Value of the velocity magnitude sensitivity. + */ + inline su2double GetTotal_Sens_ModVel() const override { return Total_Sens_ModVel; } + + /*! + * \brief Get the shape sensitivity coefficient. + * \param[in] val_marker - Surface marker where the coefficient is computed. + * \param[in] val_vertex - Vertex of the marker val_marker where the coefficient is evaluated. + * \return Value of the sensitivity coefficient. + */ + inline su2double GetCSensitivity(unsigned short val_marker, unsigned long val_vertex) const override { + return CSensitivity[val_marker][val_vertex]; + } + + /*! + * \brief Get partial derivative dIdxt. + * \param[in] iTrim - Trim variable index. + */ + su2double GetSens_dObjective_dVariables(unsigned short iTrim) const override; + + /*! + * \brief Get matrix-vector product dAdxt^T x psi. + * \param[in] iTrim - Trim variable index. + */ + su2double GetProd_dResiduals_dVariables(unsigned short iTrim) const override; + + /*! + * \brief Set the right-hand side adjoint source term. + * \param[in] iPoint - Vertex in fluid domain. + * \param[in] iVar - Variable index. + * \param[in] val - Value of the adjoint source term. + */ + void SetAdjoint_SourceTerm(unsigned long iPoint, unsigned short iVar, su2double val) override; + + /*--- Matrix-level accessors for CPyWrapperMatrixView. ---*/ + const su2matrix& GetPartialMatrix_dObjective_dStates() const override { return Partial_Sens_dObjective_dStates; } + const su2matrix& GetPartialMatrix_dResiduals_dStates() const override { return Partial_Prod_dResiduals_dStates; } + const su2matrix& GetPartialMatrix_dTractions_dStates() const override { return Partial_Prod_dTractions_dStates; } + const su2matrix& GetPartialMatrix_dObjective_dCoordinates() const override { return Partial_Sens_dObjective_dCoordinates; } + const su2matrix& GetPartialMatrix_dResiduals_dCoordinates() const override { return Partial_Prod_dResiduals_dCoordinates; } + const su2matrix& GetPartialMatrix_dTractions_dCoordinates() const override { return Partial_Prod_dTractions_dCoordinates; } + const su2matrix& GetPartialMatrix_dCoordinates_dCoordinates() const override { return Partial_Prod_dCoordinates_dCoordinates; } + const su2matrix& GetPartialMatrix_dObjective_dDisplacements(unsigned short iMarker) const override { return Partial_Sens_dObjective_dDisplacements[iMarker]; } + const su2matrix& GetPartialMatrix_dResiduals_dDisplacements(unsigned short iMarker) const override { return Partial_Prod_dResiduals_dDisplacements[iMarker]; } + const su2matrix& GetPartialMatrix_dTractions_dDisplacements(unsigned short iMarker) const override { return Partial_Prod_dTractions_dDisplacements[iMarker]; } + const su2matrix& GetPartialMatrix_dCoordinates_dDisplacements(unsigned short iMarker) const override { return Partial_Prod_dCoordinates_dDisplacements[iMarker]; } + + /*! + * \brief Prepare the solver for a new recording. + * \param[in] kind_recording - Kind of AD recording. + */ + void SetRecording(CGeometry* geometry, CConfig* config) override; + + /*! + * \brief A virtual member. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] reset - If true reset variables to their initial values. + */ + void RegisterVariables(CGeometry* geometry, CConfig* config, bool reset = false) override; + + /*! + * \brief A virtual member. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] output - Kind of output variables. + */ + void ExtractAdjoint_Variables(CGeometry* geometry, CConfig* config, ENUM_VARIABLE variable); + + /*! + * \brief A virtual member. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] output - Kind of output variables. + */ + void ExtractAdjoint_Coordinates(CGeometry* geometry, CConfig* config, CSolver* mesh_solver, ENUM_VARIABLE variable); + + /*! + * \brief A virtual member. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] config - Definition of the particular problem. + * \param[in] output - Kind of output variables. + */ + void ExtractAdjoint_Displacements(CGeometry* geometry, CConfig* config, CSolver* mesh_solver, + ENUM_VARIABLE variable); + + /*! + * \brief Update the dual-time derivatives. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver_container - Container vector with all the solutions. + * \param[in] config - Definition of the particular problem. + * \param[in] iMesh - Index of the mesh in multigrid computations. + * \param[in] iRKStep - Current step of the Runge-Kutta iteration. + * \param[in] RunTime_EqSystem - System of equations which is going to be solved. + * \param[in] Output - boolean to determine whether to print output. + */ + void Preprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config, unsigned short iMesh, + unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) override; + + /*! + * \brief Load a solution from a restart file. + * \param[in] geometry - Geometrical definition of the problem. + * \param[in] solver - Container vector with all of the solvers. + * \param[in] config - Definition of the particular problem. + * \param[in] val_iter - Current external iteration number. + * \param[in] val_update_geo - Flag for updating coords and grid velocity. + */ + void LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, + bool val_update_geo) override; + + /*! + * \brief Depends on the direct solver. + */ + inline bool GetHasHybridParallel() const override { + if (direct_solver) return direct_solver->GetHasHybridParallel(); + return false; + } +}; diff --git a/SU2_CFD/include/solvers/CSolver.hpp b/SU2_CFD/include/solvers/CSolver.hpp index eba2e5cc07d..46bbb515da1 100644 --- a/SU2_CFD/include/solvers/CSolver.hpp +++ b/SU2_CFD/include/solvers/CSolver.hpp @@ -330,6 +330,12 @@ class CSolver { */ void SetResidual_BGS(const CGeometry *geometry, const CConfig *config); + /*! + * \brief Set the value of the max residual and RMS residual. + * \param[in] val_iterlinsolver - Number of linear iterations. + */ + void ComputeResidual_RMS(const CGeometry* geometry, const CConfig* config); + /*! * \brief Set the value of the max residual and RMS residual. * \param[in] val_iterlinsolver - Number of linear iterations. @@ -1778,6 +1784,108 @@ class CSolver { CNumerics *numerics, CConfig *config) { } + /*! + * \brief Get sensitivity of deformed volume coordinates with respect to surface coordinates as a matrix-vector + * product with the adjoint variable. + * \brief Get sensitivity of objective function with respect to farfield design variables as a partial derivative. + * \param[in] iTrim - Trim variable index. + * \return Sensitivity of aerodynamic function with respect to design variable (Mach and AoA for now). + */ + inline virtual su2double GetSens_dObjective_dVariables(unsigned short iTrim) const { return su2double(0.0); } + + /*! + * \brief Get sensitivity of flow residuals with respect to farfield design variables as a matrix-vector product with + * the adjoint variable. + * \param[in] iTrim - Trim variable index. + * \return Sensitivity of aerodynamic residuals with respect to design variable (Mach and AoA for now). + */ + inline virtual su2double GetProd_dResiduals_dVariables(unsigned short iTrim) const { return su2double(0.0); } + + /*! + * \brief Set the right-hand side adjoint source term. + * \param[in] iPoint - Point index. + * \param[in] iVar - Variable index. + * \param[in] value - Value of the adjoint source term. + */ + inline virtual void SetAdjoint_SourceTerm(unsigned long iPoint, unsigned short iVar, su2double value){}; + + /*! \brief Get a reference to the matrix of partial sensitivities dObjective/dStates. */ + inline virtual const su2matrix& GetPartialMatrix_dObjective_dStates() const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + + /*! \brief Get a reference to the matrix of partial products dResiduals/dStates^T * psi. */ + inline virtual const su2matrix& GetPartialMatrix_dResiduals_dStates() const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + + /*! \brief Get a reference to the matrix of partial products dTractions/dStates^T * psi. */ + inline virtual const su2matrix& GetPartialMatrix_dTractions_dStates() const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + + /*! \brief Get a reference to the matrix of partial sensitivities dObjective/dCoordinates. */ + inline virtual const su2matrix& GetPartialMatrix_dObjective_dCoordinates() const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + + /*! \brief Get a reference to the matrix of partial products dResiduals/dCoordinates^T * psi. */ + inline virtual const su2matrix& GetPartialMatrix_dResiduals_dCoordinates() const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + + /*! \brief Get a reference to the matrix of partial products dTractions/dCoordinates^T * psi. */ + inline virtual const su2matrix& GetPartialMatrix_dTractions_dCoordinates() const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + + /*! \brief Get a reference to the matrix of partial products dCoordinates/dCoordinates^T * psi. */ + inline virtual const su2matrix& GetPartialMatrix_dCoordinates_dCoordinates() const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + + /*! \brief Get a reference to the marker matrix of partial sensitivities dObjective/dDisplacements. */ + inline virtual const su2matrix& GetPartialMatrix_dObjective_dDisplacements(unsigned short iMarker) const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + + /*! \brief Get a reference to the marker matrix of partial products dResiduals/dDisplacements^T * psi. */ + inline virtual const su2matrix& GetPartialMatrix_dResiduals_dDisplacements(unsigned short iMarker) const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + + /*! \brief Get a reference to the marker matrix of partial products dTractions/dDisplacements^T * psi. */ + inline virtual const su2matrix& GetPartialMatrix_dTractions_dDisplacements(unsigned short iMarker) const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + + /*! \brief Get a reference to the marker matrix of partial products dCoordinates/dDisplacements^T * psi. */ + inline virtual const su2matrix& GetPartialMatrix_dCoordinates_dDisplacements(unsigned short iMarker) const { + SU2_MPI::Error("Not available for this solver.", CURRENT_FUNCTION); + static su2matrix sens_matrix; + return sens_matrix; + } + /*! * \author H. Kline * \brief Provide the total "combo" objective (weighted sum of other values). @@ -3566,6 +3674,35 @@ class CSolver { */ inline virtual void ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm){} + /*! + * \brief A virtual member (overloaded for CDiscAdjResidualSolver). + * \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] output - Kind of output variables. + */ + inline virtual void ExtractAdjoint_Solution(CGeometry* geometry, CConfig* config, ENUM_VARIABLE variable) {} + + /*! + * \brief A virtual member (overloaded for CDiscAdjResidualSolver). + * \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 = false. + * \param[in] objective - Kind of output variables. + */ + inline virtual void ExtractAdjoint_Coordinates(CGeometry* geometry, CConfig* config, CSolver* mesh_solver, + ENUM_VARIABLE variable) {} + + /*! + * \brief A virtual member (overloaded for CDiscAdjResidualSolver). + * \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 = false. + * \param[in] objective - Kind of output variables. + */ + inline virtual void ExtractAdjoint_Displacements(CGeometry* geometry, CConfig* config, CSolver* mesh_solver, + ENUM_VARIABLE variable) {} + /*! * \brief Register In- or Output. * \param[in] input - Boolean whether In- or Output should be registered. @@ -3768,6 +3905,15 @@ class CSolver { */ inline virtual void ExtractAdjoint_Variables(CGeometry *geometry, CConfig *config) { } + /*! + * \brief A virtual member (overloaded for CDiscAdjResidualSolver). + * \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] output - Kind of output variables. + */ + inline virtual void ExtractAdjoint_Variables(CGeometry* geometry, CConfig* config, ENUM_VARIABLE variable) {} + /*! * \brief A virtual member. * \param[in] config - Definition of the particular problem. @@ -4281,6 +4427,21 @@ class CSolver { return VertexTraction[iMarker][iVertex][iDim]; } + /*! + * \brief Get the adjoints of the vertex tractions. + * \param[in] iMarker - Index of the marker. + * \param[in] iVertex - Index of the relevant vertex. + * \param[in] iDim - Dimension. + */ + inline su2double GetAdjointVertexTractions(unsigned short iMarker, unsigned long iVertex, unsigned short iDim) const { + return VertexTractionAdjoint[iMarker][iVertex][iDim]; + } + + /*! \brief Get a reference to the adjoint vertex traction matrix for a marker. */ + inline const su2activematrix& GetVertexTractionAdjoint(unsigned short iMarker) const { + return VertexTractionAdjoint[iMarker]; + } + /*! * \brief Register the vertex tractions as output. * \param[in] geometry - Geometrical definition. diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp index ad7255e235f..c18f097915a 100644 --- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp @@ -1,7 +1,7 @@ /*! * \file driver_adjoint_singlezone.cpp * \brief The main subroutines for driving adjoint single-zone problems. - * \author R. Sanchez + * \author R. Sanchez, H. Patel, A. Gastaldi * \version 8.4.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -32,6 +32,7 @@ #include "../../include/iteration/CIterationFactory.hpp" #include "../../include/iteration/CTurboIteration.hpp" #include "../../../Common/include/toolboxes/CQuasiNewtonInvLeastSquares.hpp" +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" CDiscAdjSinglezoneDriver::CDiscAdjSinglezoneDriver(char* confFile, unsigned short val_nZone, @@ -76,8 +77,9 @@ CDiscAdjSinglezoneDriver::CDiscAdjSinglezoneDriver(char* confFile, MainVariables = RECORDING::SOLUTION_VARIABLES; if (config->GetDeform_Mesh()) { SecondaryVariables = RECORDING::MESH_DEFORM; + } else { + SecondaryVariables = RECORDING::MESH_COORDS; } - else { SecondaryVariables = RECORDING::MESH_COORDS; } MainSolver = ADJFLOW_SOL; break; @@ -124,6 +126,10 @@ CDiscAdjSinglezoneDriver::~CDiscAdjSinglezoneDriver() { delete direct_iteration; delete direct_output; + delete PrimalJacobian; + delete PrimalPreconditioner; + delete AdjOperator; + delete AdjPreconditioner; } @@ -152,6 +158,14 @@ void CDiscAdjSinglezoneDriver::Preprocess(unsigned long TimeIter) { } void CDiscAdjSinglezoneDriver::Run() { + if (config->GetKind_DiscreteAdjoint() == ENUM_DISC_ADJ_TYPE::RESIDUALS) { + RunResidual(); + } else { + RunFixedPoint(); + } +} + +void CDiscAdjSinglezoneDriver::RunFixedPoint() { CQuasiNewtonInvLeastSquares fixPtCorrector; if (config->GetnQuasiNewtonSamples() > 1) { @@ -171,30 +185,9 @@ void CDiscAdjSinglezoneDriver::Run() { config->SetInnerIter(Adjoint_Iter); - iteration->InitializeAdjoint(solver_container, geometry_container, config_container, ZONE_0, INST_0); - - /*--- Initialize the adjoint of the objective function with 1.0. ---*/ - - SetAdjObjFunction(); - - /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ - - AD::ComputeAdjoint(); - - /*--- Extract the computed adjoint values of the input variables and store them for the next iteration. ---*/ + /*--- Update the state of the tape. ---*/ - iteration->IterateDiscAdj(geometry_container, solver_container, - config_container, ZONE_0, INST_0, false); - - /*--- Monitor the pseudo-time ---*/ - - StopCalc = iteration->Monitor(output_container[ZONE_0], integration_container, geometry_container, - solver_container, numerics_container, config_container, - surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); - - /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ - - AD::ClearAdjoints(); + UpdateAdjointsFixedPoint(); /*--- Output files for steady state simulations. ---*/ @@ -211,9 +204,208 @@ void CDiscAdjSinglezoneDriver::Run() { GetAllSolutions(ZONE_0, true, fixPtCorrector.FPresult()); SetAllSolutions(ZONE_0, true, fixPtCorrector.compute()); } + } +} + +void CDiscAdjSinglezoneDriver::RunResidual() { + if (!KrylovSet) { + /*--- Initialize the solution, right-hand-side, and system. ---*/ + const auto nVar = GetTotalNumberOfVariables(ZONE_0, true); + const auto nPoint = geometry_container[ZONE_0][INST_0][MESH_0]->GetnPoint(); + const auto nPointDomain = geometry_container[ZONE_0][INST_0][MESH_0]->GetnPointDomain(); + + AdjRHS.Initialize(nPoint, nPointDomain, nVar, nullptr); + AdjSol.Initialize(nPoint, nPointDomain, nVar, nullptr); + + AdjSolver.SetToleranceType(LinearToleranceType::RELATIVE); + KrylovSet = true; + + /*--- Initialize the preconditioner using the (transpose) approximate Jacobian from the primal problem. ---*/ + + if (config->GetKind_TimeIntScheme() != EULER_IMPLICIT) { + std::cout << "Cannot build a preconditioner for the discrete-adjoint system (missing primal Jacobian structure) !" + << std::endl; + + return; + }; + + UpdateJacobians(); + + CopiedJacobian.Initialize(nPoint, nPointDomain, nVar, nVar, true, geometry, config); + + for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { + for (unsigned long jPoint = 0; jPoint < nPoint; jPoint++) { + auto value = solver[FLOW_SOL]->Jacobian.GetBlock(iPoint, jPoint); + + CopiedJacobian.SetBlock(iPoint, jPoint, value); + } + } + + CopiedJacobian.TransposeInPlace(); + PrimalJacobian = new CSysMatrixVectorProduct(CopiedJacobian, geometry, config); + + const auto kindPreconditioner = static_cast(config->GetKind_Linear_Solver_Prec()); + PrimalPreconditioner = CPreconditioner::Create(kindPreconditioner, CopiedJacobian, geometry, config); + PrimalPreconditioner->Build(); + } + + /*--- Use FGMRES to solve the adjoint system, where: + * * the RHS is -dObjective/dStates (and any external contributions), + * * the solution are the adjoint variables, and + * * the system applies the matrix-vector product with dResidual/dStates. ---*/ + UpdateAdjointsResidual(); + + GetAllSolutions(ZONE_0, true, AdjSol); + GetAllObjectiveStatesSensitivities(AdjRHS); + + /*--- Manipulate the screen output frequency to avoid printing garbage. ---*/ + const bool monitor = true; + const auto wrtFreq = 1; + + AdjSolver.SetMonitoringFrequency(wrtFreq); + + /*--- Initialize the linear solver iterations ---*/ + AdjOperator = new LinOperator(this); + AdjPreconditioner = new LinPreconditioner(this); + + Scalar eps = 1.0; + + unsigned long nKrylov_Iter; + for (nKrylov_Iter = nAdjoint_Iter; nKrylov_Iter >= KrylovMinIters && eps > KrylovSysTol;) { + std::cout << "Adjoint iteration: " << nKrylov_Iter << " ... " << std::endl; + + auto nIter = min(nKrylov_Iter - 2ul, config->GetLinear_Solver_Iter()); + Scalar eps_l = 0.0; + Scalar tol_l = KrylovSysTol / eps; + + nIter = AdjSolver.FGMRES_LinSolver(AdjRHS, AdjSol, *AdjOperator, *AdjPreconditioner, tol_l, nIter, eps_l, monitor, + config); + nKrylov_Iter -= nIter + 1; + + eps *= eps_l; + } + + /*--- Store the solution and restore user settings. ---*/ + SetAllSolutions(ZONE_0, true, AdjSol); + + UpdateAdjointsResidual(); + + /*--- Apply the solution to obtain the total sensitivities (w.r.t. deformed volume coordinates). ---*/ + Postprocess(); + + /*--- HACK ! Force output here until proper convergence monitoring can be implemented. ---*/ + const auto inst = config_container[ZONE_0]->GetiInst(); + + for (iInst = 0; iInst < nInst[ZONE_0]; ++iInst) { + config_container[ZONE_0]->SetiInst(iInst); + output_container[ZONE_0]->SetResultFiles(geometry_container[ZONE_0][iInst][MESH_0], config_container[ZONE_0], + solver_container[ZONE_0][iInst][MESH_0], nAdjoint_Iter - nKrylov_Iter, + true); + } + config_container[ZONE_0]->SetiInst(inst); +} + +void CDiscAdjSinglezoneDriver::UpdateAdjoints() { + if (config->GetKind_DiscreteAdjoint() == ENUM_DISC_ADJ_TYPE::RESIDUALS) { + UpdateAdjointsResidual(); + } else { + UpdateAdjointsFixedPoint(); + } +} + +void CDiscAdjSinglezoneDriver::UpdateAdjointsFixedPoint() { + /*--- Initialize the adjoint of the output variables of the iteration with the adjoint solution + *--- of the previous iteration. The values are passed to the AD tool. + *--- Issues with iteration number should be dealt with once the output structure is in place. ---*/ + iteration->InitializeAdjoint(solver_container, geometry_container, config_container, ZONE_0, INST_0); + + /*--- Initialize the adjoint of the objective function with 1.0. ---*/ + SetAdjointObjective(); + + /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ + AD::ComputeAdjoint(); + + /*--- Extract the computed adjoint values of the input variables and store them for the next iteration. ---*/ + iteration->IterateDiscAdj(geometry_container, solver_container, + config_container, ZONE_0, INST_0, false); + + /*--- Monitor the pseudo-time ---*/ + StopCalc = iteration->Monitor(output_container[ZONE_0], integration_container, geometry_container, + solver_container, numerics_container, config_container, + surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); + + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ + AD::ClearAdjoints(); + +} + +void CDiscAdjSinglezoneDriver::UpdateAdjointsResidual() { + /*--- Initialize the adjoint of the output variables of the iteration with the adjoint solution + *--- of the previous iteration. The values are passed to the AD tool. + *--- Issues with iteration number should be dealt with once the output structure is in place. ---*/ + iteration->InitializeAdjoint(solver_container, geometry_container, config_container, ZONE_0, INST_0); + + /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ + AD::ComputeAdjoint(); + + /*--- Extract the adjoints of the residuals and store them for the next iteration ---*/ + if (config->GetFluidProblem()) { + solver[ADJFLOW_SOL]->ExtractAdjoint_Solution(geometry, config, ENUM_VARIABLE::RESIDUALS); + solver[ADJFLOW_SOL]->ExtractAdjoint_Variables(geometry, config, ENUM_VARIABLE::RESIDUALS); + } + + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ + AD::ClearAdjoints(); + + /*--- Initialize the adjoint of the objective function with 1.0. ---*/ + SetAdjointObjective(); + + /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ + AD::ComputeAdjoint(); + + /*--- Extract the adjoints of the objective function and store them for the next iteration ---*/ + if (config->GetFluidProblem()) { + solver[ADJFLOW_SOL]->ExtractAdjoint_Solution(geometry, config, ENUM_VARIABLE::OBJECTIVE); + solver[ADJFLOW_SOL]->ExtractAdjoint_Variables(geometry, config, ENUM_VARIABLE::OBJECTIVE); + } + + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ + AD::ClearAdjoints(); + + /*--- Initialize the adjoint of the vertex tractions with the corresponding adjoint vector. ---*/ + solver[FLOW_SOL]->SetVertexTractionsAdjoint(geometry, config); + + /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ + AD::ComputeAdjoint(); + /*--- Extract the adjoints of the vertex tractions and store them for the next iteration ---*/ + if (config->GetFluidProblem()) { + solver[ADJFLOW_SOL]->ExtractAdjoint_Solution(geometry, config, ENUM_VARIABLE::TRACTIONS); } + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ + AD::ClearAdjoints(); +} + +void CDiscAdjSinglezoneDriver::SetAdjointObjective() { + su2double seeding = 1.0; + + if (config->GetTime_Domain()) { + const auto IterAvg_Obj = config->GetIter_Avg_Objective(); + if (TimeIter < IterAvg_Obj) { + /*--- Default behavior when no window is chosen is to use Square-Windowing, i.e. the numerator equals 1.0 ---*/ + auto windowEvaluator = CWindowingTools(); + const su2double weight = windowEvaluator.GetWndWeight(config->GetKindWindow(), TimeIter, IterAvg_Obj - 1); + seeding = weight / IterAvg_Obj; + } else { + seeding = 0.0; + } + } + if (rank == MASTER_NODE) { + SU2_TYPE::SetDerivative(ObjFunc, SU2_TYPE::GetValue(seeding)); + } else { + SU2_TYPE::SetDerivative(ObjFunc, 0.0); + } } void CDiscAdjSinglezoneDriver::Postprocess() { @@ -226,6 +418,13 @@ void CDiscAdjSinglezoneDriver::Postprocess() { /*--- Compute the geometrical sensitivities ---*/ SecondaryRecording(); + + if (config->GetKind_DiscreteAdjoint() == ENUM_DISC_ADJ_TYPE::RESIDUALS) { + SecondaryRunResidual(); + } else { + SecondaryRunFixedPoint(); + } + break; case MAIN_SOLVER::DISC_ADJ_FEM : @@ -233,6 +432,12 @@ void CDiscAdjSinglezoneDriver::Postprocess() { /*--- Compute the geometrical sensitivities ---*/ SecondaryRecording(); + if (config->GetKind_DiscreteAdjoint() == ENUM_DISC_ADJ_TYPE::RESIDUALS) { + SecondaryRunResidual(); + } else { + SecondaryRunFixedPoint(); + } + iteration->Postprocess(output_container[ZONE_0], integration_container, geometry_container, solver_container, numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); @@ -288,9 +493,12 @@ void CDiscAdjSinglezoneDriver::SetRecording(RECORDING kind_recording){ iteration->SetDependencies(solver_container, geometry_container, numerics_container, config_container, ZONE_0, INST_0, kind_recording); - /*--- Do one iteration of the direct solver ---*/ - - DirectRun(kind_recording); + /*--- Do one iteration of the direct solver. ---*/ + if (config->GetKind_DiscreteAdjoint() == ENUM_DISC_ADJ_TYPE::RESIDUALS) { + DirectRunResidual(kind_recording); + } else { + DirectRunFixedPoint(kind_recording); + } /*--- Store the recording state ---*/ @@ -300,9 +508,13 @@ void CDiscAdjSinglezoneDriver::SetRecording(RECORDING kind_recording){ iteration->RegisterOutput(solver_container, geometry_container, config_container, ZONE_0, INST_0); - /*--- Extract the objective function and store it --- */ + /*--- Extract the tractions and objective function and store them. --- */ + UpdateTractions(); + UpdateObjective(); - SetObjFunction(); + if (rank == MASTER_NODE) { + AD::RegisterOutput(ObjFunc); + } if (kind_recording != RECORDING::CLEAR_INDICES && config_container[ZONE_0]->GetWrt_AD_Statistics()) { AD::PrintStatistics(SU2_MPI::GetComm(), rank == MASTER_NODE); @@ -312,72 +524,7 @@ void CDiscAdjSinglezoneDriver::SetRecording(RECORDING kind_recording){ } -void CDiscAdjSinglezoneDriver::SetAdjObjFunction(){ - su2double seeding = 1.0; - - if (config->GetTime_Domain()) { - const auto IterAvg_Obj = config->GetIter_Avg_Objective(); - if (TimeIter < IterAvg_Obj) { - /*--- Default behavior when no window is chosen is to use Square-Windowing, i.e. the numerator equals 1.0 ---*/ - auto windowEvaluator = CWindowingTools(); - const su2double weight = windowEvaluator.GetWndWeight(config->GetKindWindow(), TimeIter, IterAvg_Obj - 1); - seeding = weight / IterAvg_Obj; - } - else { - seeding = 0.0; - } - } - if (rank == MASTER_NODE) { - SU2_TYPE::SetDerivative(ObjFunc, SU2_TYPE::GetValue(seeding)); - } else { - SU2_TYPE::SetDerivative(ObjFunc, 0.0); - } -} - -void CDiscAdjSinglezoneDriver::SetObjFunction(){ - - ObjFunc = 0.0; - - /*--- Specific scalar objective functions ---*/ - - switch (config->GetKind_Solver()) { - case MAIN_SOLVER::DISC_ADJ_INC_EULER: case MAIN_SOLVER::DISC_ADJ_INC_NAVIER_STOKES: case MAIN_SOLVER::DISC_ADJ_INC_RANS: - case MAIN_SOLVER::DISC_ADJ_EULER: case MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES: case MAIN_SOLVER::DISC_ADJ_RANS: - case MAIN_SOLVER::DISC_ADJ_FEM_EULER: case MAIN_SOLVER::DISC_ADJ_FEM_NS: case MAIN_SOLVER::DISC_ADJ_FEM_RANS: - - /*--- Surface based obj. function ---*/ - - direct_output->SetHistoryOutput(geometry, solver, config, config->GetTimeIter(), - config->GetOuterIter(), config->GetInnerIter()); - ObjFunc += solver[FLOW_SOL]->GetTotal_ComboObj(); - break; - - case MAIN_SOLVER::DISC_ADJ_HEAT: - direct_output->SetHistoryOutput(geometry, solver, config, config->GetTimeIter(), - config->GetOuterIter(), config->GetInnerIter()); - ObjFunc = solver[HEAT_SOL]->GetTotal_ComboObj(); - break; - - case MAIN_SOLVER::DISC_ADJ_FEM: - solver[FEA_SOL]->Postprocessing(geometry, config, numerics_container[ZONE_0][INST_0][MESH_0][FEA_SOL], true); - - direct_output->SetHistoryOutput(geometry, solver, config, config->GetTimeIter(), - config->GetOuterIter(), config->GetInnerIter()); - ObjFunc = solver[FEA_SOL]->GetTotal_ComboObj(); - break; - - default: - break; - } - - if (rank == MASTER_NODE){ - AD::RegisterOutput(ObjFunc); - } - -} - -void CDiscAdjSinglezoneDriver::DirectRun(RECORDING kind_recording){ - +void CDiscAdjSinglezoneDriver::DirectRunFixedPoint(RECORDING kind_recording) { /*--- Mesh movement ---*/ direct_iteration->SetMesh_Deformation(geometry_container[ZONE_0][INST_0], solver, numerics, config, kind_recording); @@ -394,10 +541,28 @@ void CDiscAdjSinglezoneDriver::DirectRun(RECORDING kind_recording){ direct_iteration->Postprocess(direct_output, integration_container, geometry_container, solver_container, numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, INST_0); - /*--- Print the direct residual to screen ---*/ - + /*--- Print the direct residual to screen. ---*/ PrintDirectResidual(kind_recording); +} + +void CDiscAdjSinglezoneDriver::DirectRunResidual(RECORDING kind_recording) { + /*--- Deform the mesh. ---*/ + DeformGeometry(); + + /*--- Pre-process the primal solver state. ---*/ + UpdateTimeIter(); + UpdateFarfield(); + UpdateGeometry(); + UpdateStates(); + + /*--- Run the computation of the outputs of interest. ---*/ + UpdateResiduals(); + SetWallNormalConstraint(); + UpdateTractions(); + UpdateObjective(); + /*--- Print the direct residual to screen. ---*/ + PrintDirectResidual(kind_recording); } void CDiscAdjSinglezoneDriver::MainRecording(){ @@ -421,15 +586,16 @@ void CDiscAdjSinglezoneDriver::SecondaryRecording(){ /*--- Store the computational graph of one direct iteration with the secondary variables as input. ---*/ SetRecording(SecondaryVariables); +} +void CDiscAdjSinglezoneDriver::SecondaryRunFixedPoint() { /*--- Initialize the adjoint of the output variables of the iteration with the adjoint solution * of the current iteration. The values are passed to the AD tool. ---*/ iteration->InitializeAdjoint(solver_container, geometry_container, config_container, ZONE_0, INST_0); /*--- Initialize the adjoint of the objective function with 1.0. ---*/ - - SetAdjObjFunction(); + SetAdjointObjective(); /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ @@ -439,13 +605,271 @@ void CDiscAdjSinglezoneDriver::SecondaryRecording(){ if (SecondaryVariables == RECORDING::MESH_COORDS) { solver[MainSolver]->SetSensitivity(geometry, config); - } - else { // MESH_DEFORM + } else { solver[ADJMESH_SOL]->SetSensitivity(geometry, config, solver[MainSolver]); } + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ + AD::ClearAdjoints(); +} + +void CDiscAdjSinglezoneDriver::SecondaryRunResidual() { + /*--- Initialize the adjoint of the output variables of the iteration with the adjoint solution + *--- of the previous iteration. The values are passed to the AD tool. + *--- Issues with iteration number should be dealt with once the output structure is in place. ---*/ + + iteration->InitializeAdjoint(solver_container, geometry_container, config_container, ZONE_0, INST_0); + + /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ + + AD::ComputeAdjoint(); + + /*--- Extract the adjoints of the residuals and store them for the next iteration ---*/ + + if (config->GetFluidProblem()) { + if (SecondaryVariables == RECORDING::MESH_COORDS) { + solver[ADJFLOW_SOL]->ExtractAdjoint_Coordinates(geometry, config, nullptr, ENUM_VARIABLE::RESIDUALS); + } else { + solver[ADJFLOW_SOL]->ExtractAdjoint_Coordinates(geometry, config, solver[ADJMESH_SOL], ENUM_VARIABLE::RESIDUALS); + } + } + + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ + + AD::ClearAdjoints(); + + /*--- Initialize the adjoint of the objective function with 1.0. ---*/ + + SetAdjointObjective(); + + /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ + + AD::ComputeAdjoint(); + + /*--- Extract the adjoints of the objective function and store them for the next iteration ---*/ + + if (config->GetFluidProblem()) { + if (SecondaryVariables == RECORDING::MESH_COORDS) { + solver[ADJFLOW_SOL]->ExtractAdjoint_Coordinates(geometry, config, nullptr, ENUM_VARIABLE::OBJECTIVE); + } else { + solver[ADJFLOW_SOL]->ExtractAdjoint_Coordinates(geometry, config, solver[ADJMESH_SOL], ENUM_VARIABLE::OBJECTIVE); + } + } + + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ + + AD::ClearAdjoints(); + + /*--- Initialize the adjoint of the vertex tractions with the corresponding adjoint vector. ---*/ + + solver[FLOW_SOL]->SetVertexTractionsAdjoint(geometry, config); + + /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ + + AD::ComputeAdjoint(); + + /*--- Extract the adjoints of the vertex tractions and store them for the next iteration ---*/ + + if (config->GetFluidProblem()) { + if (SecondaryVariables == RECORDING::MESH_COORDS) { + solver[ADJFLOW_SOL]->ExtractAdjoint_Coordinates(geometry, config, nullptr, ENUM_VARIABLE::TRACTIONS); + } else { + solver[ADJFLOW_SOL]->ExtractAdjoint_Coordinates(geometry, config, solver[ADJMESH_SOL], ENUM_VARIABLE::TRACTIONS); + } + } + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ AD::ClearAdjoints(); + /*--- Skip the derivation of the mesh solver if it is not defined ---*/ + if (SecondaryVariables == RECORDING::MESH_DEFORM) { + /*--- Initialize the adjoint of the volume coordinates with the corresponding adjoint vector. ---*/ + SU2_OMP_PARALLEL_(if(solver[ADJMESH_SOL]->GetHasHybridParallel())) { + + /*--- Initialize the adjoints of the volume coordinates ---*/ + solver[ADJMESH_SOL]->SetAdjoint_Output(geometry, config); + } + END_SU2_OMP_PARALLEL + + /*--- Interpret the stored information by calling the corresponding routine of the AD tool. ---*/ + AD::ComputeAdjoint(); + + /*--- Extract the adjoints of the volume coordinates and store them for the next iteration ---*/ + if (config->GetFluidProblem()) { + solver[ADJFLOW_SOL]->ExtractAdjoint_Coordinates(geometry, config, solver[ADJMESH_SOL], + ENUM_VARIABLE::COORDINATES); + } + + /*--- Clear the stored adjoint information to be ready for a new evaluation. ---*/ + AD::ClearAdjoints(); + } + + /*--- Extract the adjoints of the residuals and store them for the next iteration ---*/ + if (config->GetFluidProblem()) { + solver[ADJFLOW_SOL]->SetSensitivity(geometry, config); + } +} + +void CDiscAdjSinglezoneDriver::UpdateTimeIter() { + /*--- Update the primal time iteration. ---*/ + solver[FLOW_SOL]->SetTime_Step(geometry, solver, config, MESH_0, config->GetTimeIter()); +} + +void CDiscAdjSinglezoneDriver::UpdateFarfield() { + /*--- Update the primal far-field variables. ---*/ + + su2double Velocity_Ref = config->GetVelocity_Ref(); + su2double Alpha = config->GetAoA() * PI_NUMBER / 180.0; + su2double Beta = config->GetAoS() * PI_NUMBER / 180.0; + su2double Mach = config->GetMach(); + su2double Temperature = config->GetTemperature_FreeStream(); + su2double Gas_Constant = config->GetGas_Constant(); + su2double Gamma = config->GetGamma(); + su2double SoundSpeed = sqrt(Gamma * Gas_Constant * Temperature); + + if (nDim == 2) { + config->GetVelocity_FreeStreamND()[0] = cos(Alpha) * Mach * SoundSpeed / Velocity_Ref; + config->GetVelocity_FreeStreamND()[1] = sin(Alpha) * Mach * SoundSpeed / Velocity_Ref; + } + if (nDim == 3) { + config->GetVelocity_FreeStreamND()[0] = cos(Alpha) * cos(Beta) * Mach * SoundSpeed / Velocity_Ref; + config->GetVelocity_FreeStreamND()[1] = sin(Beta) * Mach * SoundSpeed / Velocity_Ref; + config->GetVelocity_FreeStreamND()[2] = sin(Alpha) * cos(Beta) * Mach * SoundSpeed / Velocity_Ref; + } +} + +void CDiscAdjSinglezoneDriver::UpdateGeometry() { + /*--- Update the geometry (i.e. dual grid, without multi-grid). ---*/ + + geometry->InitiateComms(geometry, config, MPI_QUANTITIES::COORDINATES); + geometry->CompleteComms(geometry, config, MPI_QUANTITIES::COORDINATES); + + geometry->SetControlVolume(config, UPDATE); + geometry->SetBoundControlVolume(config, UPDATE); + geometry->SetMaxLength(config); +} + +void CDiscAdjSinglezoneDriver::DeformGeometry() { + /*--- Deform the geometry. ---*/ + + direct_iteration->SetMesh_Deformation(geometry_container[ZONE_0][INST_0], solver, numerics, config, + SecondaryVariables); +} + +void CDiscAdjSinglezoneDriver::UpdateObjective() { + ObjFunc = 0.0; + + /*--- Specific scalar objective functions ---*/ + switch (config->GetKind_Solver()) { + case MAIN_SOLVER::DISC_ADJ_INC_EULER: case MAIN_SOLVER::DISC_ADJ_INC_NAVIER_STOKES: case MAIN_SOLVER::DISC_ADJ_INC_RANS: + case MAIN_SOLVER::DISC_ADJ_EULER: case MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES: case MAIN_SOLVER::DISC_ADJ_RANS: + case MAIN_SOLVER::DISC_ADJ_FEM_EULER: case MAIN_SOLVER::DISC_ADJ_FEM_NS: case MAIN_SOLVER::DISC_ADJ_FEM_RANS: + + /*--- Surface based obj. function ---*/ + direct_output->SetHistoryOutput(geometry, solver, config, config->GetTimeIter(), + config->GetOuterIter(), config->GetInnerIter()); + ObjFunc += solver[FLOW_SOL]->GetTotal_ComboObj(); + break; + + case MAIN_SOLVER::DISC_ADJ_HEAT: + direct_output->SetHistoryOutput(geometry, solver, config, config->GetTimeIter(), + config->GetOuterIter(), config->GetInnerIter()); + ObjFunc = solver[HEAT_SOL]->GetTotal_ComboObj(); + break; + + case MAIN_SOLVER::DISC_ADJ_FEM: + solver[FEA_SOL]->Postprocessing(geometry, config, numerics_container[ZONE_0][INST_0][MESH_0][FEA_SOL], true); + + direct_output->SetHistoryOutput(geometry, solver, config, config->GetTimeIter(), + config->GetOuterIter(), config->GetInnerIter()); + ObjFunc = solver[FEA_SOL]->GetTotal_ComboObj(); + break; + + default: + break; + } +} + +void CDiscAdjSinglezoneDriver::UpdateStates() { + /*--- Update the flow and turbulent conservative state variables, preparing for other updates. ---*/ + direct_iteration->Preprocess(direct_output, integration_container, geometry_container, solver_container, + numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, + INST_0); +} + +void CDiscAdjSinglezoneDriver::UpdateResiduals() { + /*--- Update the Euler, Navier-Stokes or Reynolds-averaged Navier-Stokes (RANS) residuals and objective function (no + * system solve). ---*/ + integration[FLOW_SOL]->ComputeResiduals(geometry_container, solver_container, numerics_container, config_container, + FLOW_SOL, ZONE_0, INST_0); +} + +void CDiscAdjSinglezoneDriver::UpdateTractions() { + /*--- Update the surface tractions. ---*/ + direct_iteration->Postprocess(direct_output, integration_container, geometry_container, solver_container, + numerics_container, config_container, surface_movement, grid_movement, FFDBox, ZONE_0, + INST_0); +} + +void CDiscAdjSinglezoneDriver::UpdateJacobians() { + /*--- Compute the approximate Jacobian for preconditioning. ---*/ + solver[FLOW_SOL]->PrepareImplicitIteration(geometry, solver, config); +} + +void CDiscAdjSinglezoneDriver::ApplyPreconditioner(const CSysVector& u, CSysVector& v) { + /*--- Use an approximate diagonal preconditioning based on the transpose of the primal Jacobian. ---*/ + (*PrimalPreconditioner)(u, v); + + /*--- Apply a few FGMRES iterations in addition to the above preconditioner. ---*/ + Scalar KrylovPreEps = KrylovPreTol; + auto MaxIter = 5; + solver[FLOW_SOL]->System.FGMRES_LinSolver(u, v, *PrimalJacobian, *PrimalPreconditioner, KrylovPreTol, MaxIter, + KrylovPreEps, false, config); +} + +void CDiscAdjSinglezoneDriver::ApplyOperator(const CSysVector& u, CSysVector& v) { + /*--- Set the adjoint variables used in the seeding of the tape. ---*/ + SetAllSolutions(ZONE_0, true, u); + + /*--- Evaluate the tape to and extract the partial derivatives. ---*/ + UpdateAdjointsResidual(); + + /*--- Extract the partial residual Jacobian-adjoint product. ---*/ + GetAllResidualsStatesSensitivities(v); +} + +void CDiscAdjSinglezoneDriver::SetWallNormalConstraint() { + const unsigned short iVel = 1; + + for (auto iMarker = 0ul; iMarker < config->GetnMarker_All(); iMarker++) { + if (config->GetMarker_All_KindBC(iMarker) != EULER_WALL && + config->GetMarker_All_KindBC(iMarker) != SYMMETRY_PLANE) continue; + + for (auto iVertex = 0ul; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + if (!geometry->nodes->GetDomain(iPoint)) continue; + + /*--- Get the unit normal. ---*/ + const auto Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + su2double UnitNormal[3] = {0.0}; + const auto it = geometry->symmetryNormals[iMarker].find(iVertex); + + if (it != geometry->symmetryNormals[iMarker].end()) { + for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] = it->second[iDim]; + } else { + const su2double Area = GeometryToolbox::Norm(nDim, Normal); + for (auto iDim = 0u; iDim < nDim; iDim++) UnitNormal[iDim] = Normal[iDim] / Area; + } + + /*--- Compute wall-normal momentum from the registered (unmodified) solution. ---*/ + su2double vnMom = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) + vnMom += solver[FLOW_SOL]->GetNodes()->GetSolution(iPoint, iVel + iDim) * UnitNormal[iDim]; + + /*--- Add the constraint R_n = (ρv)·n into the currently zeroed normal slot. ---*/ + for (auto iDim = 0u; iDim < nDim; iDim++) + solver[FLOW_SOL]->LinSysRes(iPoint, iVel + iDim) += vnMom * UnitNormal[iDim]; + } + } } diff --git a/SU2_CFD/src/integration/CIntegration.cpp b/SU2_CFD/src/integration/CIntegration.cpp index edb4248d8c1..b506cbbf069 100644 --- a/SU2_CFD/src/integration/CIntegration.cpp +++ b/SU2_CFD/src/integration/CIntegration.cpp @@ -265,3 +265,65 @@ void CIntegration::SetDualTime_Solver(const CGeometry *geometry, CSolver *solver } END_SU2_OMP_PARALLEL } + +void CIntegration::ComputeResiduals(CGeometry**** geometry_container, CSolver***** solvers_container, + CNumerics****** numerics_container, CConfig** config_container, + unsigned short EqSystem, unsigned short iZone, unsigned short iInst) { + CGeometry* geometry = geometry_container[iZone][iInst][MESH_0]; + CConfig* config = config_container[iZone]; + CSolver** solvers = solvers_container[iZone][iInst][MESH_0]; + CNumerics** numerics = numerics_container[iZone][iInst][MESH_0][EqSystem]; + + /*--- Update global parameters. ---*/ + + switch (config->GetKind_Solver()) { + case MAIN_SOLVER::EULER: + case MAIN_SOLVER::DISC_ADJ_EULER: + case MAIN_SOLVER::INC_EULER: + case MAIN_SOLVER::DISC_ADJ_INC_EULER: + case MAIN_SOLVER::NEMO_EULER: + config->SetGlobalParam(MAIN_SOLVER::EULER, RUNTIME_FLOW_SYS); + break; + + case MAIN_SOLVER::NAVIER_STOKES: + case MAIN_SOLVER::DISC_ADJ_NAVIER_STOKES: + case MAIN_SOLVER::INC_NAVIER_STOKES: + case MAIN_SOLVER::DISC_ADJ_INC_NAVIER_STOKES: + case MAIN_SOLVER::NEMO_NAVIER_STOKES: + config->SetGlobalParam(MAIN_SOLVER::NAVIER_STOKES, RUNTIME_FLOW_SYS); + break; + + case MAIN_SOLVER::RANS: + case MAIN_SOLVER::DISC_ADJ_RANS: + case MAIN_SOLVER::INC_RANS: + case MAIN_SOLVER::DISC_ADJ_INC_RANS: + config->SetGlobalParam(MAIN_SOLVER::RANS, RUNTIME_FLOW_SYS); + break; + + default: + break; + } + + /*--- Pre-processing re-sets the residuals before each evaluation. ---*/ + + solvers[EqSystem]->Preprocessing(geometry, solvers, config, MESH_0, 0, EqSystem, false); + + /*--- Space integration computes the residuals. ---*/ + + /*--- Set the Jacobian to zero since this is not done inside the fluid iteration + when running the discrete adjoint solver. ---*/ + + solvers[EqSystem]->Jacobian.SetValZero(); // TODO: Check if actually necessary (already done in Preprocessing) + + Space_Integration(geometry, solvers, numerics, config, MESH_0, NO_RK_ITER, EqSystem); + + solvers[EqSystem]->ComputeResidual_RMS(geometry, config); + + /*--- Post-processing. ---*/ + + solvers[EqSystem]->Postprocessing(geometry, solvers, config, MESH_0); + + solvers[EqSystem]->Pressure_Forces(geometry, config); + solvers[EqSystem]->Momentum_Forces(geometry, config); + solvers[EqSystem]->Friction_Forces(geometry, config); +} diff --git a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp index a10f7c4505f..371cc418948 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFluidIteration.cpp @@ -540,7 +540,9 @@ void CDiscAdjFluidIteration::RegisterOutput(CSolver***** solver, CGeometry**** g if (config[iZone]->GetFluidProblem() && !config[iZone]->GetMultizone_Problem()) { solvers0[FLOW_SOL]->RegisterVertexTractions(geometry0, config[iZone]); } - + if (config[iZone]->GetDeform_Mesh()) { + solver[iZone][iInst][MESH_0][ADJMESH_SOL]->RegisterOutput(geometry[iZone][iInst][MESH_0], config[iZone]); + } } END_SU2_OMP_PARALLEL } diff --git a/SU2_CFD/src/meson.build b/SU2_CFD/src/meson.build index 3b40822a34f..10975b6443d 100644 --- a/SU2_CFD/src/meson.build +++ b/SU2_CFD/src/meson.build @@ -97,6 +97,7 @@ su2_cfd_src += files(['solvers/CSolverFactory.cpp', 'solvers/CRadSolver.cpp', 'solvers/CRadP1Solver.cpp', 'solvers/CDiscAdjSolver.cpp', + 'solvers/CDiscAdjResidualSolver.cpp', 'solvers/CEulerSolver.cpp', 'solvers/CFEASolver.cpp', 'solvers/CFEM_DG_EulerSolver.cpp', diff --git a/SU2_CFD/src/python_wrapper_structure.cpp b/SU2_CFD/src/python_wrapper_structure.cpp index 8c5f939305d..c3366658517 100644 --- a/SU2_CFD/src/python_wrapper_structure.cpp +++ b/SU2_CFD/src/python_wrapper_structure.cpp @@ -53,7 +53,206 @@ void CDriver::PreprocessPythonInterface(CConfig** config, CGeometry**** geometry } ////////////////////////////////////////////////////////////////////////////////// -/* Functions to obtain global parameters from SU2 (time steps, delta t, etc.) */ +/* Functions related to farfield flow variables. */ +////////////////////////////////////////////////////////////////////////////////// +passivedouble CDriver::GetAngleOfAttack() const { return SU2_TYPE::GetValue(main_config->GetAoA()); } + +void CDriver::SetAngleOfAttack(const passivedouble AoA) { + config_container[selected_zone]->SetAoA(AoA); + solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL]->UpdateFarfieldVelocity(config_container[selected_zone]); +} + +passivedouble CDriver::GetAngleOfSideslip() const { return SU2_TYPE::GetValue(main_config->GetAoS()); } + +void CDriver::SetAngleOfSideslip(const passivedouble AoS) { + config_container[selected_zone]->SetAoS(AoS); + solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL]->UpdateFarfieldVelocity(config_container[selected_zone]); +} + +passivedouble CDriver::GetMachNumber() const { return SU2_TYPE::GetValue(main_config->GetMach()); } + +void CDriver::SetMachNumber(passivedouble value) { + main_config->SetMach(value); + UpdateFarfield(); +} + +passivedouble CDriver::GetReynoldsNumber() const { return SU2_TYPE::GetValue(main_config->GetReynolds()); } + +void CDriver::SetReynoldsNumber(passivedouble value) { + main_config->SetReynolds(value); + UpdateFarfield(); +} + +////////////////////////////////////////////////////////////////////////////////// +/* Functions related to the flow solver solution and variables. */ +////////////////////////////////////////////////////////////////////////////////// + +unsigned long CDriver::GetNumberStateVariables() const { + if (!main_config->GetFluidProblem()) { + SU2_MPI::Error("Flow solver is not defined!", CURRENT_FUNCTION); + } + + return solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetnVar(); +} + +unsigned long CDriver::GetNumberPrimitiveVariables() const { + if (!main_config->GetFluidProblem()) { + SU2_MPI::Error("Flow solver is not defined!", CURRENT_FUNCTION); + } + + return solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetnPrimVar(); +} + +////////////////////////////////////////////////////////////////////////////////// +/* Functions related to the adjoint flow solver solution. */ +////////////////////////////////////////////////////////////////////////////////// + +CPyWrapperMatrixView CDriver::MarkerAdjointForces(unsigned short iMarker) const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast(solver_container[ZONE_0][INST_0][MESH_0][FLOW_SOL]->GetVertexTractionAdjoint(iMarker)); + return CPyWrapperMatrixView(mat, "MarkerAdjointForces", true); +} + +CPyWrapperMatrixView CDriver::CoordinatesCoordinatesSensitivities() const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dCoordinates_dCoordinates()); + return CPyWrapperMatrixView(mat, "CoordinatesCoordinatesSensitivities", true); +} + +CPyWrapperMatrixView CDriver::MarkerCoordinatesDisplacementsSensitivities(unsigned short iMarker) const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dCoordinates_dDisplacements(iMarker)); + return CPyWrapperMatrixView(mat, "MarkerCoordinatesDisplacementsSensitivities", true); +} + +vector CDriver::GetObjectiveFarfieldVariablesSensitivities() const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + if (main_config->GetKind_DiscreteAdjoint() != ENUM_DISC_ADJ_TYPE::RESIDUALS) { + SU2_MPI::Error("Discrete adjoint flow solver does not use residual-based formulation!", CURRENT_FUNCTION); + } + + const int nTrim = 2; + vector values(nTrim, 0.0); + + su2double mach = solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetSens_dObjective_dVariables(0); + values[0] = SU2_TYPE::GetValue(mach); + + su2double alpha = solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetSens_dObjective_dVariables(1); + values[1] = SU2_TYPE::GetValue(alpha); + + return values; +} + +vector CDriver::GetResidualsFarfieldVariablesSensitivities() const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + if (main_config->GetKind_DiscreteAdjoint() != ENUM_DISC_ADJ_TYPE::RESIDUALS) { + SU2_MPI::Error("Discrete adjoint flow solver does not use residual-based formulation!", CURRENT_FUNCTION); + } + + const int nTrim = 2; + vector values(nTrim, 0.0); + + su2double mach = solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetProd_dResiduals_dVariables(0); + values[0] = SU2_TYPE::GetValue(mach); + + su2double alpha = solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetProd_dResiduals_dVariables(1); + values[1] = SU2_TYPE::GetValue(alpha); + + return values; +} + +CPyWrapperMatrixView CDriver::ObjectiveStatesSensitivities() const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dObjective_dStates()); + return CPyWrapperMatrixView(mat, "ObjectiveStatesSensitivities", true); +} + +CPyWrapperMatrixView CDriver::ResidualsStatesSensitivities() const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dResiduals_dStates()); + return CPyWrapperMatrixView(mat, "ResidualsStatesSensitivities", true); +} + +CPyWrapperMatrixView CDriver::ForcesStatesSensitivities() const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dTractions_dStates()); + return CPyWrapperMatrixView(mat, "ForcesStatesSensitivities", true); +} + +CPyWrapperMatrixView CDriver::ObjectiveCoordinatesSensitivities() const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dObjective_dCoordinates()); + return CPyWrapperMatrixView(mat, "ObjectiveCoordinatesSensitivities", true); +} + +CPyWrapperMatrixView CDriver::ResidualsCoordinatesSensitivities() const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dResiduals_dCoordinates()); + return CPyWrapperMatrixView(mat, "ResidualsCoordinatesSensitivities", true); +} + +CPyWrapperMatrixView CDriver::ForcesCoordinatesSensitivities() const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dTractions_dCoordinates()); + return CPyWrapperMatrixView(mat, "ForcesCoordinatesSensitivities", true); +} + +CPyWrapperMatrixView CDriver::MarkerObjectiveDisplacementsSensitivities(unsigned short iMarker) const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dObjective_dDisplacements(iMarker)); + return CPyWrapperMatrixView(mat, "MarkerObjectiveDisplacementsSensitivities", true); +} + +CPyWrapperMatrixView CDriver::MarkerResidualsDisplacementsSensitivities(unsigned short iMarker) const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dResiduals_dDisplacements(iMarker)); + return CPyWrapperMatrixView(mat, "MarkerResidualsDisplacementsSensitivities", true); +} + +CPyWrapperMatrixView CDriver::MarkerForcesDisplacementsSensitivities(unsigned short iMarker) const { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dTractions_dDisplacements(iMarker)); + return CPyWrapperMatrixView(mat, "MarkerForcesDisplacementsSensitivities", true); +} + +CPyWrapperMatrixView CDriver::AdjointSourceTerm() { + if (!main_config->GetFluidProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint flow solver is not defined!", CURRENT_FUNCTION); + } + auto& mat = const_cast&>(solver_container[ZONE_0][INST_0][MESH_0][ADJFLOW_SOL]->GetPartialMatrix_dObjective_dStates()); + return CPyWrapperMatrixView(mat, "AdjointSourceTerm", false); +} + +////////////////////////////////////////////////////////////////////////////////// +/* Functions to obtain global parameters from SU2 (time steps, delta t, etc.). */ ////////////////////////////////////////////////////////////////////////////////// unsigned long CDriver::GetNumberTimeIter() const { return config_container[selected_zone]->GetnTime_Iter(); } @@ -74,9 +273,22 @@ passivedouble CDriver::GetUnsteadyTimeStep() const { string CDriver::GetSurfaceFileName() const { return config_container[selected_zone]->GetSurfCoeff_FileName(); } -//////////////////////////////////////////////////////////////////////////////// -/* Functions related to the management of markers */ -//////////////////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////////////////////// +/* Functions related to the management of markers. */ +////////////////////////////////////////////////////////////////////////////////// + +vector CDriver::GetFluidLoadMarkerTags() const { + const auto nMarker = main_config->GetnMarker_Fluid_Load(); + vector tags; + + tags.resize(nMarker); + + for (auto iMarker = 0u; iMarker < nMarker; iMarker++) { + tags[iMarker] = main_config->GetMarker_Fluid_Load_TagBound(iMarker); + } + + return tags; +} void CDriver::SetHeatSourcePosition(passivedouble alpha, passivedouble pos_x, passivedouble pos_y, passivedouble pos_z) { @@ -101,19 +313,9 @@ void CDriver::SetInletAngle(unsigned short iMarker, passivedouble alpha) { } } -void CDriver::SetFarFieldAoA(const passivedouble AoA) { - config_container[selected_zone]->SetAoA(AoA); - solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL]->UpdateFarfieldVelocity(config_container[selected_zone]); -} - -void CDriver::SetFarFieldAoS(const passivedouble AoS) { - config_container[selected_zone]->SetAoS(AoS); - solver_container[selected_zone][INST_0][MESH_0][FLOW_SOL]->UpdateFarfieldVelocity(config_container[selected_zone]); -} - -///////////////////////////////////////////////////////////////////////////////////////////////////////////////// -/* Functions related to simulation control, high level functions (reset convergence, set initial mesh, etc.) */ -///////////////////////////////////////////////////////////////////////////////////////////////////////////////// +//////////////////////////////////////////////////////////////////////////////////////////////////////////////// +/* Functions related to simulation control, high level functions (reset convergence, set initial mesh, etc.). */ +//////////////////////////////////////////////////////////////////////////////////////////////////////////////// void CSinglezoneDriver::SetInitialMesh() { DynamicMeshUpdate(0); @@ -151,9 +353,66 @@ void CDriver::BoundaryConditionsUpdate() { } } -//////////////////////////////////////////////////////////////////////////////// -/* Functions related to dynamic mesh */ -//////////////////////////////////////////////////////////////////////////////// +void CDriver::UpdateGeometry() { + geometry_container[ZONE_0][INST_0][MESH_0]->InitiateComms(main_geometry, main_config, MPI_QUANTITIES::COORDINATES); + geometry_container[ZONE_0][INST_0][MESH_0]->CompleteComms(main_geometry, main_config, MPI_QUANTITIES::COORDINATES); + + geometry_container[ZONE_0][INST_0][MESH_0]->SetControlVolume(main_config, UPDATE); + geometry_container[ZONE_0][INST_0][MESH_0]->SetBoundControlVolume(main_config, UPDATE); + geometry_container[ZONE_0][INST_0][MESH_0]->SetMaxLength(main_config); +} + +void CDriver::UpdateFarfield() { + su2double Velocity_Ref = main_config->GetVelocity_Ref(); + su2double Alpha = main_config->GetAoA() * PI_NUMBER / 180.0; + su2double Beta = main_config->GetAoS() * PI_NUMBER / 180.0; + su2double Mach = main_config->GetMach(); + su2double Temperature = main_config->GetTemperature_FreeStream(); + su2double Gas_Constant = main_config->GetGas_Constant(); + su2double Gamma = main_config->GetGamma(); + su2double SoundSpeed = sqrt(Gamma * Gas_Constant * Temperature); + + if (nDim == 2) { + main_config->GetVelocity_FreeStreamND()[0] = cos(Alpha) * Mach * SoundSpeed / Velocity_Ref; + main_config->GetVelocity_FreeStreamND()[1] = sin(Alpha) * Mach * SoundSpeed / Velocity_Ref; + } + if (nDim == 3) { + main_config->GetVelocity_FreeStreamND()[0] = cos(Alpha) * cos(Beta) * Mach * SoundSpeed / Velocity_Ref; + main_config->GetVelocity_FreeStreamND()[1] = sin(Beta) * Mach * SoundSpeed / Velocity_Ref; + main_config->GetVelocity_FreeStreamND()[2] = sin(Alpha) * Mach * SoundSpeed / Velocity_Ref; + } +} + +////////////////////////////////////////////////////////////////////////////////// +/* Functions related to adjoint finite element simulations. */ +////////////////////////////////////////////////////////////////////////////////// + +vector CDriver::GetMarkerForceSensitivities(unsigned short iMarker) const { + if (!main_config->GetStructuralProblem() || !main_config->GetDiscrete_Adjoint()) { + SU2_MPI::Error("Discrete adjoint structural solver is not defined!", CURRENT_FUNCTION); + } + if (main_config->GetKind_DiscreteAdjoint() != ENUM_DISC_ADJ_TYPE::FIXED_POINT) { + SU2_MPI::Error("Discrete adjoint structural solver does not use fixed-point formulation!", CURRENT_FUNCTION); + } + + const auto nVertex = GetNumberMarkerNodes(iMarker); + vector values(nVertex * nDim, 0.0); + + for (auto iVertex = 0ul; iVertex < nVertex; iVertex++) { + auto iPoint = main_geometry->vertex[iMarker][iVertex]->GetNode(); + + for (auto iDim = 0u; iDim < nDim; iDim++) { + values[iVertex * nDim + iDim] = SU2_TYPE::GetValue( + solver_container[ZONE_0][INST_0][MESH_0][ADJFEA_SOL]->GetNodes()->GetFlowTractionSensitivity(iPoint, iDim)); + } + } + + return values; +} + +////////////////////////////////////////////////////////////////////////////////// +/* Functions related to dynamic mesh. */ +////////////////////////////////////////////////////////////////////////////////// void CDriver::SetTranslationRate(passivedouble xDot, passivedouble yDot, passivedouble zDot) { main_config->SetTranslation_Rate(0, xDot); @@ -178,4 +437,3 @@ void CDriver::SetMarkerTranslationRate(unsigned short iMarker, passivedouble vel config_container[selected_zone]->SetMarkerTranslationRate(iMarker, 1, vel_y); config_container[selected_zone]->SetMarkerTranslationRate(iMarker, 2, vel_z); } - diff --git a/SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp index a7c37f0ece0..d03183b2f0e 100644 --- a/SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CDiscAdjMeshSolver.cpp @@ -97,11 +97,16 @@ void CDiscAdjMeshSolver::SetRecording(CGeometry* geometry, CConfig *config){ } +void CDiscAdjMeshSolver::RegisterOutput(CGeometry *geometry, CConfig *config) { + + /*--- Register variables as output of the solver iteration.---*/ + direct_solver->GetNodes()->RegisterSolution(false); +} + void CDiscAdjMeshSolver::RegisterSolution(CGeometry *geometry, CConfig *config){ /*--- Register reference mesh coordinates ---*/ direct_solver->GetNodes()->Register_MeshCoord(); - } void CDiscAdjMeshSolver::RegisterVariables(CGeometry *geometry, CConfig *config, bool reset){ @@ -114,6 +119,19 @@ void CDiscAdjMeshSolver::RegisterVariables(CGeometry *geometry, CConfig *config, SU2_OMP_SAFE_GLOBAL_ACCESS(direct_solver->GetNodes()->Register_BoundDisp();) } +void CDiscAdjMeshSolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config){ + + su2double Solution[MAXNVAR] = {0.0}; + unsigned short iVar; + + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++){ + for (iVar = 0; iVar < nVar; iVar++) + Solution[iVar] = nodes->GetSolution(iPoint,iVar); + + direct_solver->GetNodes()->SetAdjointSolution(iPoint,Solution); + } +} + void CDiscAdjMeshSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm){ /*--- Extract the sensitivities of the mesh coordinates ---*/ @@ -147,7 +165,7 @@ void CDiscAdjMeshSolver::ExtractAdjoint_Variables(CGeometry *geometry, CConfig * su2double Solution[MAXNVAR] = {0.0}; direct_solver->GetNodes()->GetAdjoint_BoundDisp(iPoint,Solution); - nodes->SetBoundDisp_Sens(iPoint,Solution); + nodes->SetBoundDisp_Sens(iPoint,Solution); } END_SU2_OMP_FOR diff --git a/SU2_CFD/src/solvers/CDiscAdjResidualSolver.cpp b/SU2_CFD/src/solvers/CDiscAdjResidualSolver.cpp new file mode 100755 index 00000000000..553ed595a1f --- /dev/null +++ b/SU2_CFD/src/solvers/CDiscAdjResidualSolver.cpp @@ -0,0 +1,693 @@ +/*! + * \file CDiscAdjResidualSolver.cpp + * \brief Main subroutines for solving the residual-based discrete adjoint problem. + * \author H. Patel, A. Gastaldi + * \version 8.1.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../include/solvers/CDiscAdjResidualSolver.hpp" + +#include "../../../Common/include/toolboxes/geometry_toolbox.hpp" + +CDiscAdjResidualSolver::CDiscAdjResidualSolver(CGeometry *geometry, CConfig *config, CSolver *direct_solver, + unsigned short Kind_Solver, unsigned short iMesh) : CSolver() { + adjoint = true; + + nVar = direct_solver->GetnVar(); + nDim = geometry->GetnDim(); + + /*--- Initialize arrays to NULL. ---*/ + + /*-- Store some information about direct solver. ---*/ + this->KindDirect_Solver = Kind_Solver; + this->direct_solver = direct_solver; + + nMarker = config->GetnMarker_All(); + nPoint = geometry->GetnPoint(); + nPointDomain = geometry->GetnPointDomain(); + + omp_chunk_size = computeStaticChunkSize(nPoint, omp_get_max_threads(), OMP_MAX_SIZE); + + /*--- Define some auxiliary vectors related to the residual. ---*/ + + Residual_RMS.resize(nVar, 1.0); + Residual_Max.resize(nVar, 1.0); + Point_Max.resize(nVar, 0); + Point_Max_Coord.resize(nVar, nDim) = su2double(0.0); + + /*--- Define some auxiliary vectors related to the residual for problems with a BGS strategy. ---*/ + + if (config->GetMultizone_Residual()) { + Residual_BGS.resize(nVar, 1.0); + Residual_Max_BGS.resize(nVar, 1.0); + Point_Max_BGS.resize(nVar, 0); + Point_Max_Coord_BGS.resize(nVar, nDim) = su2double(0.0); + } + + /*--- Define some auxiliary vectors related to the residual for problems with a BGS strategy. ---*/ + + // Derivatives w.r.t. flow variables + Partial_Prod_dResiduals_dVariables.resize(nTrim) = su2double(0.0); + Partial_Sens_dObjective_dVariables.resize(nTrim) = su2double(0.0); + + // Derivatives w.r.t. flow states + Partial_Prod_dResiduals_dStates.resize(nPoint, nVar) = su2double(0.0); + Partial_Sens_dObjective_dStates.resize(nPoint, nVar) = su2double(0.0); + Partial_Prod_dTractions_dStates.resize(nPoint, nVar) = su2double(0.0); + + // Derivatives w.r.t. mesh coordinates + Partial_Prod_dResiduals_dCoordinates.resize(nPoint, nDim) = su2double(0.0); + Partial_Sens_dObjective_dCoordinates.resize(nPoint, nDim) = su2double(0.0); + Partial_Prod_dTractions_dCoordinates.resize(nPoint, nDim) = su2double(0.0); + + Partial_Prod_dCoordinates_dCoordinates.resize(nPoint, nDim) = su2double(0.0); + + // Derivatives w.r.t. mesh displacements on each marker + Partial_Prod_dResiduals_dDisplacements.resize(nMarker); + Partial_Sens_dObjective_dDisplacements.resize(nMarker); + Partial_Prod_dTractions_dDisplacements.resize(nMarker); + + Partial_Prod_dCoordinates_dDisplacements.resize(nMarker); + + for (auto iMarker = 0ul; iMarker < nMarker; iMarker++) { + const auto nVertex = geometry->GetnVertex(iMarker); + + Partial_Prod_dResiduals_dDisplacements[iMarker].resize(nVertex, nDim) = su2double(0.0); + Partial_Sens_dObjective_dDisplacements[iMarker].resize(nVertex, nDim) = su2double(0.0); + Partial_Prod_dTractions_dDisplacements[iMarker].resize(nVertex, nDim) = su2double(0.0); + + Partial_Prod_dCoordinates_dDisplacements[iMarker].resize(nVertex, nDim) = su2double(0.0); + } + + AD_ResidualIndex.resize(nPoint, nVar) = AD::Identifier(-1); + + /*--- Sensitivity definition and coefficient in all the markers. ---*/ + + CSensitivity.resize(nMarker); + for (auto iMarker = 0ul; iMarker < nMarker; iMarker++) { + const auto nVertex = geometry->nVertex[iMarker]; + CSensitivity[iMarker].resize(nVertex, 0.0); + } + + Sens_Geo.resize(config->GetnMarker_Monitoring(), 0.0); + + /*--- Initialize the discrete adjoint solution to zero everywhere. ---*/ + + if (nVar > MAXNVAR) { + SU2_MPI::Error("Oops! The CDiscAdjResidualSolver static array sizes are not large enough.", CURRENT_FUNCTION); + } + + vector Solution(nVar, 1e-16); + nodes = new CDiscAdjVariable(Solution.data(), nPoint, nDim, nVar, config); + SetBaseClassPointerToNodes(); + + switch (KindDirect_Solver) { + case RUNTIME_FLOW_SYS: + SolverName = "ADJ.FLOW"; + break; + case RUNTIME_HEAT_SYS: + SolverName = "ADJ.HEAT"; + break; + case RUNTIME_TURB_SYS: + SolverName = "ADJ.TURB"; + break; + case RUNTIME_SPECIES_SYS: + SolverName = "ADJ.SPECIES"; + break; + case RUNTIME_RADIATION_SYS: + SolverName = "ADJ.RAD"; + break; + default: + SolverName = "ADJ.SOL"; + break; + } +} + +CDiscAdjResidualSolver::~CDiscAdjResidualSolver(void) { delete nodes; } + +void CDiscAdjResidualSolver::SetRecording(CGeometry* geometry, CConfig* config) { + const bool time_n1_needed = config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND; + const bool time_n_needed = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || time_n1_needed; + + /*--- Reset the solution to the initial (converged) solution. ---*/ + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + direct_solver->GetNodes()->SetSolution(iPoint, nodes->GetSolution_Direct(iPoint)); + } + END_SU2_OMP_FOR + + if (time_n_needed) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + for (auto iVar = 0u; iVar < nVar; iVar++) { + AD::ResetInput(direct_solver->GetNodes()->GetSolution_time_n(iPoint)[iVar]); + } + } + END_SU2_OMP_FOR + } + if (time_n1_needed) { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + for (auto iVar = 0u; iVar < nVar; iVar++) { + AD::ResetInput(direct_solver->GetNodes()->GetSolution_time_n1(iPoint)[iVar]); + } + } + END_SU2_OMP_FOR + } + + /*--- Set indices to zero. ---*/ + + RegisterVariables(geometry, config, true); +} + +void CDiscAdjResidualSolver::RegisterSolution(CGeometry* geometry, CConfig* config) { + const bool time_n1_needed = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); + const bool time_n_needed = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || time_n1_needed; + + /*--- Register solution at all necessary time instances and other variables on the tape. ---*/ + + /*--- Boolean true indicates that an input is registered. ---*/ + direct_solver->GetNodes()->RegisterSolution(true); + + if (time_n_needed) direct_solver->GetNodes()->RegisterSolution_time_n(); + + if (time_n1_needed) direct_solver->GetNodes()->RegisterSolution_time_n1(); +} + +void CDiscAdjResidualSolver::RegisterVariables(CGeometry* geometry, CConfig* config, bool reset) { + SU2_OMP_MASTER { + /*--- Register farfield values as input. ---*/ + + if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && + (KindDirect_Solver == RUNTIME_FLOW_SYS && !config->GetBoolTurbomachinery())) { + su2double Velocity_Ref = config->GetVelocity_Ref(); + Alpha = config->GetAoA() * PI_NUMBER / 180.0; + Beta = config->GetAoS() * PI_NUMBER / 180.0; + Mach = config->GetMach(); + Pressure = config->GetPressure_FreeStreamND(); + Temperature = config->GetTemperature_FreeStreamND(); + + su2double SoundSpeed = 0.0; + + if (nDim == 2) { + SoundSpeed = config->GetVelocity_FreeStreamND()[0] * Velocity_Ref / (cos(Alpha) * Mach); + } + if (nDim == 3) { + SoundSpeed = config->GetVelocity_FreeStreamND()[0] * Velocity_Ref / (cos(Alpha) * cos(Beta) * Mach); + } + + if (!reset) { + AD::RegisterInput(Mach); + AD::RegisterInput(Alpha); + AD::RegisterInput(Temperature); + AD::RegisterInput(Pressure); + } + + /*--- Recompute the free-stream velocity. ---*/ + + if (nDim == 2) { + config->GetVelocity_FreeStreamND()[0] = cos(Alpha) * Mach * SoundSpeed / Velocity_Ref; + config->GetVelocity_FreeStreamND()[1] = sin(Alpha) * Mach * SoundSpeed / Velocity_Ref; + } + if (nDim == 3) { + config->GetVelocity_FreeStreamND()[0] = cos(Alpha) * cos(Beta) * Mach * SoundSpeed / Velocity_Ref; + config->GetVelocity_FreeStreamND()[1] = sin(Beta) * Mach * SoundSpeed / Velocity_Ref; + config->GetVelocity_FreeStreamND()[2] = sin(Alpha) * cos(Beta) * Mach * SoundSpeed / Velocity_Ref; + } + // TODO: The following is not necessary if modifying Mach and Alpha of solver, rather than config + // OR: Alternatively, the config values could be registered instead of the solver values + config->SetAoA(Alpha * 180.0 / PI_NUMBER); + config->SetMach(Mach); + + config->SetTemperature_FreeStreamND(Temperature); + direct_solver->SetTemperature_Inf(Temperature); + config->SetPressure_FreeStreamND(Pressure); + direct_solver->SetPressure_Inf(Pressure); + } + + if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS) && + config->GetBoolTurbomachinery()) { + BPressure = config->GetPressureOut_BC(); + Temperature = config->GetTotalTemperatureIn_BC(); + + if (!reset) { + AD::RegisterInput(BPressure); + AD::RegisterInput(Temperature); + } + + config->SetPressureOut_BC(BPressure); + config->SetTotalTemperatureIn_BC(Temperature); + } + + /*--- Register incompressible initialization values as input. ---*/ + + if ((config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) && + ((KindDirect_Solver == RUNTIME_FLOW_SYS && (!config->GetBoolTurbomachinery())))) { + /*--- Access the velocity (or pressure) and temperature at the + inlet BC and the back pressure at the outlet. Note that we are + assuming that have internal flow, which will be true for the + majority of cases. External flows with far-field BCs will report + zero for these sensitivities. ---*/ + + ModVel = config->GetIncInlet_BC(); + BPressure = config->GetIncPressureOut_BC(); + Temperature = config->GetIncTemperature_BC(); + + /*--- Register the variables for AD. ---*/ + + if (!reset) { + AD::RegisterInput(ModVel); + AD::RegisterInput(BPressure); + AD::RegisterInput(Temperature); + } + + /*--- Set the BC values in the config class. ---*/ + + config->SetIncInlet_BC(ModVel); + config->SetIncPressureOut_BC(BPressure); + config->SetIncTemperature_BC(Temperature); + } + + /*--- Register incompressible radiation values as input. ---*/ + + if ((config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE) && + ((KindDirect_Solver == RUNTIME_RADIATION_SYS && (!config->GetBoolTurbomachinery())))) { + + /*--- Access the nondimensional free-stream temperature. ---*/ + + TemperatureRad = config->GetTemperature_FreeStreamND(); + + /*--- Register the variables for AD. ---*/ + + if (!reset) { + AD::RegisterInput(TemperatureRad); + } + + /*--- Set the temperature at infinity in the direct solver class. ---*/ + + direct_solver->SetTemperature_Inf(TemperatureRad); + } + + /*--- Here it is possible to register other variables as input that influence the flow solution + and thereby also the objective function. The adjoint values (i.e. the derivatives) can be + extracted in the ExtractAdjointVariables routine. ---*/ + } + END_SU2_OMP_MASTER + SU2_OMP_BARRIER +} + +void CDiscAdjResidualSolver::RegisterOutput(CGeometry* geometry, CConfig* config) { + /*--- Register variables as output of the solver iteration. Boolean false indicates that an output is registered. ---*/ + + SU2_OMP_FOR_STAT(roundUpDiv(nPoint, omp_get_num_threads())) + for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) { + for (unsigned long iVar = 0; iVar < nVar; ++iVar) { + AD::RegisterOutput(direct_solver->LinSysRes(iPoint, iVar)); + + auto& residual_index = AD_ResidualIndex(iPoint, iVar); + AD::SetIndex(residual_index, direct_solver->LinSysRes(iPoint, iVar)); + } + } + END_SU2_OMP_FOR +} + +void CDiscAdjResidualSolver::ExtractAdjoint_Solution(CGeometry* geometry, CConfig* config, ENUM_VARIABLE variable) { + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + if (variable == ENUM_VARIABLE::OBJECTIVE) { + direct_solver->GetNodes()->GetAdjointSolution(iPoint, Partial_Sens_dObjective_dStates[iPoint]); + } else if (variable == ENUM_VARIABLE::RESIDUALS) { + direct_solver->GetNodes()->GetAdjointSolution(iPoint, Partial_Prod_dResiduals_dStates[iPoint]); + } else if (variable == ENUM_VARIABLE::TRACTIONS) { + direct_solver->GetNodes()->GetAdjointSolution(iPoint, Partial_Prod_dTractions_dStates[iPoint]); + } else { + SU2_MPI::Error("The discrete adjoint solver does not support this as an output variable.", CURRENT_FUNCTION); + } + } +} + +void CDiscAdjResidualSolver::ExtractAdjoint_Variables(CGeometry* geometry, CConfig* config, ENUM_VARIABLE variable) { + SU2_OMP_MASTER { + /*--- Extract the adjoint values of the farfield values. ---*/ + + if ((config->GetKind_Regime() == ENUM_REGIME::COMPRESSIBLE) && (KindDirect_Solver == RUNTIME_FLOW_SYS) && + !config->GetBoolTurbomachinery()) { + su2double Local_Sens_AoA, Local_Sens_Mach; + + Local_Sens_Mach = SU2_TYPE::GetDerivative(Mach); + Local_Sens_AoA = SU2_TYPE::GetDerivative(Alpha); + + SU2_MPI::Allreduce(&Local_Sens_Mach, &Total_Sens_Mach, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + SU2_MPI::Allreduce(&Local_Sens_AoA, &Total_Sens_AoA, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + } else { + SU2_MPI::Error("The discrete adjoint solver does not support this as an output variable.", CURRENT_FUNCTION); + } + + if (variable == ENUM_VARIABLE::OBJECTIVE) { + Partial_Sens_dObjective_dVariables[0] = Total_Sens_Mach; + Partial_Sens_dObjective_dVariables[1] = Total_Sens_AoA; + } else if (variable == ENUM_VARIABLE::RESIDUALS) { + Partial_Prod_dResiduals_dVariables[0] = Total_Sens_Mach; + Partial_Prod_dResiduals_dVariables[1] = Total_Sens_AoA; + } else { + SU2_MPI::Error("The discrete adjoint solver does not support this as an output variable.", CURRENT_FUNCTION); + } + } + END_SU2_OMP_MASTER + SU2_OMP_BARRIER +} + +void CDiscAdjResidualSolver::ExtractAdjoint_Coordinates(CGeometry* geometry, CConfig* config, CSolver* mesh_solver, + ENUM_VARIABLE variable) { + SU2_OMP_PARALLEL { + su2matrix* VolumeSensitivity = nullptr; + + if (variable == ENUM_VARIABLE::OBJECTIVE) { + VolumeSensitivity = &Partial_Sens_dObjective_dCoordinates; + } else if (variable == ENUM_VARIABLE::RESIDUALS) { + VolumeSensitivity = &Partial_Prod_dResiduals_dCoordinates; + } else if (variable == ENUM_VARIABLE::TRACTIONS) { + VolumeSensitivity = &Partial_Prod_dTractions_dCoordinates; + } else { + SU2_MPI::Error("The discrete adjoint solver does not support this as an output variable.", CURRENT_FUNCTION); + } + + /*--- Sensitivities w.r.t all reference mesh coordinates. ---*/ + + if (mesh_solver) { + mesh_solver->ExtractAdjoint_Solution(geometry, config, true); + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + for (auto iDim = 0u; iDim < nDim; iDim++) { + (*VolumeSensitivity)(iPoint, iDim) = mesh_solver->GetNodes()->GetSolution(iPoint, iDim); + } + } + END_SU2_OMP_FOR + } else { + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + auto coord = geometry->nodes->GetCoord(iPoint); + + for (auto iDim = 0u; iDim < nDim; iDim++) { + (*VolumeSensitivity)(iPoint, iDim) = geometry->nodes->GetAdjointSolution(iPoint, iDim); + AD::ResetInput(coord[iDim]); + } + } + END_SU2_OMP_FOR + } + } + END_SU2_OMP_PARALLEL +} + +void CDiscAdjResidualSolver::ExtractAdjoint_Displacements(CGeometry* geometry, CConfig* config, CSolver* mesh_solver, + ENUM_VARIABLE variable) { + SU2_OMP_PARALLEL { + vector>* MarkerSensitivity = nullptr; + + if (variable == ENUM_VARIABLE::COORDINATES) { + MarkerSensitivity = &Partial_Prod_dCoordinates_dDisplacements; + } else if (variable == ENUM_VARIABLE::OBJECTIVE) { + MarkerSensitivity = &Partial_Sens_dObjective_dDisplacements; + } else if (variable == ENUM_VARIABLE::RESIDUALS) { + MarkerSensitivity = &Partial_Prod_dResiduals_dDisplacements; + } else if (variable == ENUM_VARIABLE::TRACTIONS) { + MarkerSensitivity = &Partial_Prod_dTractions_dDisplacements; + } else { + SU2_MPI::Error("The discrete adjoint solver does not support this as an output variable.", CURRENT_FUNCTION); + } + + /*--- If the mesh solver is used, extract the discrete-adjoint sensitivities of the boundary displacements. ---*/ + + const auto eps = config->GetAdjSharp_LimiterCoeff() * config->GetRefElemLength(); + + if (mesh_solver) { + mesh_solver->ExtractAdjoint_Variables(geometry, config); + + for (auto iMarker = 0ul; iMarker < nMarker; iMarker++) { + if (!config->GetSolid_Wall(iMarker)) continue; + + /*--- Sensitivities w.r.t aerodynamic boundary displacements. ---*/ + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0ul; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + + /*--- If sharp edge, set the sensitivity to 0. ---*/ + + su2double limiter; + if (config->GetSens_Remove_Sharp() && (geometry->nodes->GetSharpEdge_Distance(iPoint) < eps)) { + limiter = 0.0; + } else { + limiter = 1.0; + } + + /*--- Get the gradient from the deformed coordinates in the mesh solver. ---*/ + + for (auto iDim = 0u; iDim < nDim; iDim++) { + (*MarkerSensitivity)[iMarker](iVertex, iDim) = + limiter * mesh_solver->GetNodes()->GetBoundDisp_Sens(iPoint, iDim); + } + } + END_SU2_OMP_FOR + } + } + + else { + // TODO: Implement derivation of the legacy mesh deformation (see SU2_DEF_AD) + } + } + END_SU2_OMP_PARALLEL +} + +void CDiscAdjResidualSolver::SetAdjoint_Output(CGeometry* geometry, CConfig* config) { + const bool dual_time = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST || + config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); + const bool multizone = config->GetMultizone_Problem(); + + /*--- Local container to manipulate the adjoint solution. ---*/ + su2double Solution[MAXNVAR] = {0.0}; + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + /*--- Get and store the adjoint solution of a point. ---*/ + + for (auto iVar = 0u; iVar < nVar; iVar++) { + Solution[iVar] = nodes->GetSolution(iPoint, iVar); + } + + /*--- Add dual time contributions to the adjoint solution. Two terms stored for DT-2nd-order. ---*/ + + if (dual_time && !multizone) { + for (auto iVar = 0u; iVar < nVar; iVar++) { + Solution[iVar] += nodes->GetDual_Time_Derivative(iPoint, iVar); + } + } + + /*--- Set the adjoint values of the primal solution. ---*/ + + for (unsigned long iVar = 0; iVar < nVar; iVar++) + AD::SetDerivative(AD_ResidualIndex(iPoint, iVar), SU2_TYPE::GetValue(Solution[iVar])); + } + END_SU2_OMP_FOR +} + +void CDiscAdjResidualSolver::SetSensitivity(CGeometry* geometry, CConfig* config, CSolver*) { + SU2_OMP_PARALLEL { + const bool time_stepping = (config->GetTime_Marching() != TIME_MARCHING::STEADY); + const su2double eps = config->GetAdjSharp_LimiterCoeff() * config->GetRefElemLength(); + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) { + for (auto iDim = 0u; iDim < nDim; iDim++) { + su2double Sensitivity = + Partial_Sens_dObjective_dCoordinates(iPoint, iDim) + Partial_Prod_dResiduals_dCoordinates(iPoint, iDim); + + /*--- If sharp edge, set the sensitivity to 0 on that region. ---*/ + + if (config->GetSens_Remove_Sharp() && geometry->nodes->GetSharpEdge_Distance(iPoint) < eps) { + Sensitivity = 0.0; + } + + if (!time_stepping) { + nodes->SetSensitivity(iPoint, iDim, Sensitivity); + } else { + nodes->SetSensitivity(iPoint, iDim, nodes->GetSensitivity(iPoint, iDim) + Sensitivity); + } + } + } + END_SU2_OMP_FOR + + SetSurface_Sensitivity(geometry, config); + } + END_SU2_OMP_PARALLEL +} + +void CDiscAdjResidualSolver::SetSurface_Sensitivity(CGeometry* geometry, CConfig* config) { + SU2_OMP_MASTER + for (auto& x : Sens_Geo) x = 0.0; + END_SU2_OMP_MASTER + + /*--- Loop over boundary markers to select those for Euler walls and Navier-Stokes walls. ---*/ + + for (auto iMarker = 0ul; iMarker < nMarker; iMarker++) { + if (!config->GetSolid_Wall(iMarker)) continue; + + su2double Sens = 0.0; + + SU2_OMP_FOR_STAT(OMP_MIN_SIZE) + for (auto iVertex = 0ul; iVertex < geometry->GetnVertex(iMarker); iVertex++) { + /*--- Projection of the gradient calculated with AD onto the normal vector of the surface. ---*/ + + const auto iPoint = geometry->vertex[iMarker][iVertex]->GetNode(); + const auto Normal = geometry->vertex[iMarker][iVertex]->GetNormal(); + + su2double Sens_Vertex = 0.0; + for (auto iDim = 0u; iDim < nDim; iDim++) { + Sens_Vertex += Normal[iDim] * nodes->GetSensitivity(iPoint, iDim); + } + Sens_Vertex /= GeometryToolbox::Norm(nDim, Normal); + + CSensitivity[iMarker][iVertex] = -Sens_Vertex; + Sens += pow(Sens_Vertex, 2); + } + END_SU2_OMP_FOR + + if (config->GetMarker_All_Monitoring(iMarker) == NO) continue; + + /*--- Compute sensitivity for each surface point. ---*/ + + const auto Marker_Tag = config->GetMarker_All_TagBound(iMarker); + + for (size_t iMarker_Mon = 0; iMarker_Mon < Sens_Geo.size(); iMarker_Mon++) { + if (Marker_Tag == config->GetMarker_Monitoring_TagBound(iMarker_Mon)) { + atomicAdd(Sens, Sens_Geo[iMarker_Mon]); + break; + } + } + } + + SU2_OMP_BARRIER + SU2_OMP_MASTER { + auto local = Sens_Geo; + SU2_MPI::Allreduce(local.data(), Sens_Geo.data(), Sens_Geo.size(), MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm()); + + Total_Sens_Geo = 0.0; + for (auto& x : Sens_Geo) { + x = sqrt(x); + Total_Sens_Geo += x; + } + } + END_SU2_OMP_MASTER + SU2_OMP_BARRIER +} + +su2double CDiscAdjResidualSolver::GetSens_dObjective_dVariables(unsigned short iTrim) const { + return Partial_Sens_dObjective_dVariables(iTrim); +} + +su2double CDiscAdjResidualSolver::GetProd_dResiduals_dVariables(unsigned short iTrim) const { + return Partial_Prod_dResiduals_dVariables(iTrim); +} + +void CDiscAdjResidualSolver::SetAdjoint_SourceTerm(unsigned long iPoint, unsigned short iVar, su2double val) { + Partial_Sens_dObjective_dStates[iPoint][iVar] = val; +} + +void CDiscAdjResidualSolver::Preprocessing(CGeometry* geometry, CSolver** solver_container, CConfig* config, + unsigned short iMesh, unsigned short iRKStep, + unsigned short RunTime_EqSystem, bool Output) { + config->SetGlobalParam(config->GetKind_Solver(), RunTime_EqSystem); + + const bool dual_time = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) || + (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); + + if (!dual_time) return; + + SU2_OMP_FOR_STAT(omp_chunk_size) + for (auto iPoint = 0ul; iPoint < geometry->GetnPoint(); iPoint++) { + const auto solution_n = nodes->GetSolution_time_n(iPoint); + const auto solution_n1 = nodes->GetSolution_time_n1(iPoint); + + for (auto iVar = 0u; iVar < nVar; iVar++) { + nodes->SetDual_Time_Derivative(iPoint, iVar, solution_n[iVar] + nodes->GetDual_Time_Derivative_n(iPoint, iVar)); + nodes->SetDual_Time_Derivative_n(iPoint, iVar, solution_n1[iVar]); + } + END_SU2_OMP_FOR + } +} + +void CDiscAdjResidualSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig* config, int val_iter, + bool val_update_geo) { + /*--- Restart the solution from file information. ---*/ + + auto filename = config->GetSolution_AdjFileName(); + auto restart_filename = config->GetObjFunc_Extension(filename); + restart_filename = config->GetFilename(restart_filename, "", val_iter); + + const bool rans = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); + + /*--- Skip coordinates. ---*/ + + unsigned short skipVars = geometry[MESH_0]->GetnDim(); + + /*--- Skip flow adjoint variables. ---*/ + + if (KindDirect_Solver == RUNTIME_TURB_SYS) { + skipVars += nDim + 2; + } + + if (KindDirect_Solver == RUNTIME_SPECIES_SYS) { + // Skip the number of flow variables and turbulent variables to get to the adjoint species variables. + skipVars += nDim + 2; + if (rans) skipVars += solver[MESH_0][TURB_SOL]->GetnVar(); + } + + /*--- Skip flow adjoint and turbulent variables. ---*/ + + if (KindDirect_Solver == RUNTIME_RADIATION_SYS) { + skipVars += nDim + 2; + if (rans) skipVars += solver[MESH_0][TURB_SOL]->GetnVar(); + } + + BasicLoadRestart(geometry[MESH_0], config, restart_filename, skipVars); + + /*--- Interpolate solution on coarse grids. ---*/ + + for (auto iMesh = 1u; iMesh <= config->GetnMGLevels(); iMesh++) { + const auto& fineSol = solver[iMesh - 1][ADJFLOW_SOL]->GetNodes()->GetSolution(); + + for (auto iPoint = 0ul; iPoint < geometry[iMesh]->GetnPoint(); iPoint++) { + su2double Solution[MAXNVAR] = {0.0}; + const su2double Area_Parent = geometry[iMesh]->nodes->GetVolume(iPoint); + + for (auto iChildren = 0u; iChildren < geometry[iMesh]->nodes->GetnChildren_CV(iPoint); iChildren++) { + const auto Point_Fine = geometry[iMesh]->nodes->GetChildren_CV(iPoint, iChildren); + const su2double weight = geometry[iMesh - 1]->nodes->GetVolume(Point_Fine) / Area_Parent; + + for (auto iVar = 0u; iVar < nVar; iVar++) Solution[iVar] += weight * fineSol(Point_Fine, iVar); + } + solver[iMesh][ADJFLOW_SOL]->GetNodes()->SetSolution(iPoint, Solution); + } + } +} diff --git a/SU2_CFD/src/solvers/CSolver.cpp b/SU2_CFD/src/solvers/CSolver.cpp index 58f8c90e2bc..d834ddbde73 100644 --- a/SU2_CFD/src/solvers/CSolver.cpp +++ b/SU2_CFD/src/solvers/CSolver.cpp @@ -4164,6 +4164,52 @@ void CSolver::SetVerificationSolution(unsigned short nDim, } } +void CSolver::ComputeResidual_RMS(const CGeometry *geometry, const CConfig *config){ + + /*--- Set Residuals to zero ---*/ + SU2_OMP_MASTER + for (unsigned short iVar = 0; iVar < nVar; iVar++){ + Residual_RMS[iVar] = 0.0; + Residual_Max[iVar] = 0.0; + } + END_SU2_OMP_MASTER + + vector resRMS(nVar, 0.0); + vector resMax(nVar,0.0); + vector coordMax(nVar,nullptr); + vector idxMax(nVar,0); + + /*--- Iterate over residuals and determine RMS and MAX values ---*/ + + SU2_OMP_FOR_(schedule(static, OMP_MIN_SIZE) SU2_NOWAIT) + for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) { + + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + unsigned long total_index = iPoint*nVar + iVar; + + su2double Res = fabs(LinSysRes[total_index]); + resRMS[iVar] += Res*Res; + + if (Res > resMax[iVar]) { + resMax[iVar] = Res; + idxMax[iVar] = iPoint; + coordMax[iVar] = geometry->nodes->GetCoord(iPoint); + } + } + } + END_SU2_OMP_FOR + SU2_OMP_CRITICAL + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + Residual_RMS[iVar] += resRMS[iVar]; + AddRes_Max(iVar, resMax[iVar], geometry->nodes->GetGlobalIndex(idxMax[iVar]), coordMax[iVar]); + } + END_SU2_OMP_CRITICAL + SU2_OMP_BARRIER + + /*--- Compute the root mean square residual ---*/ + SetResidual_RMS(geometry, config); +} + void CSolver::ComputeResidual_Multizone(const CGeometry *geometry, const CConfig *config){ SU2_OMP_PARALLEL { diff --git a/SU2_CFD/src/solvers/CSolverFactory.cpp b/SU2_CFD/src/solvers/CSolverFactory.cpp index 31a31495367..83af09ca158 100644 --- a/SU2_CFD/src/solvers/CSolverFactory.cpp +++ b/SU2_CFD/src/solvers/CSolverFactory.cpp @@ -43,6 +43,7 @@ #include "../../include/solvers/CFEASolver.hpp" #include "../../include/solvers/CTemplateSolver.hpp" #include "../../include/solvers/CDiscAdjSolver.hpp" +#include "../../include/solvers/CDiscAdjResidualSolver.hpp" #include "../../include/solvers/CDiscAdjFEASolver.hpp" #include "../../include/solvers/CFEM_DG_EulerSolver.hpp" #include "../../include/solvers/CFEM_DG_NSSolver.hpp" @@ -261,7 +262,12 @@ CSolver* CSolverFactory::CreateSubSolver(SUB_SOLVER_TYPE kindSolver, CSolver **s metaData.integrationType = INTEGRATION_TYPE::DEFAULT; break; case SUB_SOLVER_TYPE::DISC_ADJ_FLOW: - genericSolver = new CDiscAdjSolver(geometry, config, solver[FLOW_SOL], RUNTIME_FLOW_SYS, iMGLevel); + if (config->GetKind_DiscreteAdjoint() == ENUM_DISC_ADJ_TYPE::RESIDUALS) { + genericSolver = new CDiscAdjResidualSolver(geometry, config, solver[FLOW_SOL], RUNTIME_FLOW_SYS, iMGLevel); + } + else { + genericSolver = new CDiscAdjSolver(geometry, config, solver[FLOW_SOL], RUNTIME_FLOW_SYS, iMGLevel); + } metaData.integrationType = INTEGRATION_TYPE::DEFAULT; break; case SUB_SOLVER_TYPE::EULER: diff --git a/SU2_CFD/src/variables/CDiscAdjMeshBoundVariable.cpp b/SU2_CFD/src/variables/CDiscAdjMeshBoundVariable.cpp index d84f3fc949b..74e20e70f7f 100644 --- a/SU2_CFD/src/variables/CDiscAdjMeshBoundVariable.cpp +++ b/SU2_CFD/src/variables/CDiscAdjMeshBoundVariable.cpp @@ -29,14 +29,7 @@ #include "../../include/variables/CDiscAdjMeshBoundVariable.hpp" -CDiscAdjMeshBoundVariable::CDiscAdjMeshBoundVariable(unsigned long npoint, unsigned long ndim, CConfig *config) { - - nPoint = npoint; - nVar = ndim; - nDim = ndim; - - /*--- Allocate the solution array. ---*/ - Solution.resize(nPoint,nVar) = su2double(0.0); +CDiscAdjMeshBoundVariable::CDiscAdjMeshBoundVariable(unsigned long npoint, unsigned long ndim, CConfig *config): CVariable(npoint, ndim, ndim, config, true) { VertexMap.Reset(nPoint); } diff --git a/SU2_CFD/src/variables/CMeshBoundVariable.cpp b/SU2_CFD/src/variables/CMeshBoundVariable.cpp index b906eed3665..d26f0b00fcd 100644 --- a/SU2_CFD/src/variables/CMeshBoundVariable.cpp +++ b/SU2_CFD/src/variables/CMeshBoundVariable.cpp @@ -31,6 +31,12 @@ CMeshBoundVariable::CMeshBoundVariable(unsigned long npoint, unsigned long ndim, CMeshVariable(npoint, ndim, config) { VertexMap.Reset(nPoint); + + /*--- Allocate the AD indices ---*/ + if (config->GetDiscrete_Adjoint()) { + AD_InputIndex.resize(nPoint,nDim) = AD::GetPassiveIndex(); + AD_OutputIndex.resize(nPoint,nDim) = AD::GetPassiveIndex(); + } } void CMeshBoundVariable::AllocateBoundaryVariables(CConfig *config) { diff --git a/SU2_DEF/include/drivers/CDiscAdjDeformationDriver.hpp b/SU2_DEF/include/drivers/CDiscAdjDeformationDriver.hpp index d590ca255b2..0f4e99c1f41 100644 --- a/SU2_DEF/include/drivers/CDiscAdjDeformationDriver.hpp +++ b/SU2_DEF/include/drivers/CDiscAdjDeformationDriver.hpp @@ -59,6 +59,34 @@ class CDiscAdjDeformationDriver : public CDriverBase { */ void Finalize() override; + /*! + * \brief Get total sensitivity of the objective function w.r.t. surface coordinates or displacements + * at mesh vertices. + * \return A CPyWrapperMatrixView of shape (nPoint, nDim) with the coordinate sensitivities. + */ + CPyWrapperMatrixView GetObjectiveCoordinatesTotalSensitivities() const; + + /*! + * \brief Get total sensitivity of the objective function w.r.t. surface coordinates or displacements + * at marker vertices. + * \param[in] iMarker - Marker identifier. + * \return A CPyWrapperMarkerMatrixView of shape (nVertex, nDim) with the coordinate sensitivities. + */ + CPyWrapperMarkerMatrixView GetMarkerObjectiveCoordinatesTotalSensitivities(unsigned short iMarker) const; + + /*! + * \brief Get total sensitivity of the objective function w.r.t. the design variables. + * \return Total sensitivity of the objective function w.r.t. the design variables. + */ + vector> GetObjectiveDesignVariablesTotalSensitivities() const; + + /*! + * \brief Get total sensitivity of the objective function w.r.t. a design variable. + * \param[in] iDV - Design Variable index. + * \return Total sensitivity of the objective function w.r.t. a design variable. + */ + vector GetObjectiveDesignVariablesTotalSensitivities(unsigned short iDV) const; + protected: /*! * \brief Read in the config and mesh files. diff --git a/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp b/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp index c609c956f41..0645ed17f5d 100644 --- a/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp +++ b/SU2_DEF/src/drivers/CDiscAdjDeformationDriver.cpp @@ -52,7 +52,7 @@ CDiscAdjDeformationDriver::CDiscAdjDeformationDriver(char* confFile, SU2_Comm MP Gradient = new su2double*[nDV]; for (auto iDV = 0u; iDV < nDV; iDV++) { - /*--- Initialize to zero ---*/ + /*--- structure to store the gradient. ---*/ unsigned short nValue = config_container[ZONE_0]->GetnDV_Value(iDV); Gradient[iDV] = new su2double[nValue]; @@ -290,10 +290,11 @@ void CDiscAdjDeformationDriver::Preprocess() { /*--- Read in sensitivities from file. ---*/ - if (config_container[ZONE_0]->GetSensitivity_Format() == UNORDERED_ASCII) + if (config_container[ZONE_0]->GetSensitivity_Format() == UNORDERED_ASCII) { geometry_container[iZone][INST_0][MESH_0]->ReadUnorderedSensitivity(config_container[iZone]); - else + } else { geometry_container[iZone][INST_0][MESH_0]->SetSensitivity(config_container[iZone]); + } } } } @@ -303,6 +304,7 @@ void CDiscAdjDeformationDriver::Run() { if (!config_container[iZone]->GetDiscrete_Adjoint()) { continue; } + if (rank == MASTER_NODE) cout << "\n---------------------- Mesh sensitivity computation ---------------------" << endl; if (config_container[iZone]->GetDiscrete_Adjoint() && config_container[iZone]->GetSmoothGradient() && @@ -379,7 +381,7 @@ void CDiscAdjDeformationDriver::Finalize() { /*--- Exit the solver cleanly. ---*/ if (rank == MASTER_NODE) - cout << "\n------------------------- Exit Success (SU2_DOT) ------------------------" << endl << endl; + cout << "\n------------------------- Exit Success (SU2_DEF_AD) ------------------------" << endl << endl; } void CDiscAdjDeformationDriver::SetProjection_FD(CGeometry* geometry, CConfig* config, @@ -931,6 +933,7 @@ void CDiscAdjDeformationDriver::DerivativeTreatment_MeshSensitivity(CGeometry* g if (rank == MASTER_NODE) cout << " working on mesh level" << endl; /*--- Work with the surface derivatives. ---*/ + if (config->GetSmoothOnSurface()) { /*--- Project to surface and then smooth on surface. ---*/ @@ -1001,3 +1004,44 @@ void CDiscAdjDeformationDriver::DerivativeTreatment_Gradient(CGeometry* geometry solver->RecordTapeAndCalculateOriginalGradient(geometry, surface_movement, grid_movement, config, Gradient); } } + +CPyWrapperMatrixView CDiscAdjDeformationDriver::GetObjectiveCoordinatesTotalSensitivities() const { + auto& sensitivity = const_cast(geometry_container[ZONE_0][INST_0][MESH_0]->GetSensitivityMatrix()); + return CPyWrapperMatrixView(sensitivity, "ObjectiveCoordinatesTotalSensitivities", true); +} + +CPyWrapperMarkerMatrixView CDiscAdjDeformationDriver::GetMarkerObjectiveCoordinatesTotalSensitivities( + unsigned short iMarker) const { + if (iMarker >= GetNumberMarkers()) SU2_MPI::Error("Marker index exceeds size.", CURRENT_FUNCTION); + auto& sensitivity = const_cast(geometry_container[ZONE_0][INST_0][MESH_0]->GetSensitivityMatrix()); + return CPyWrapperMarkerMatrixView(sensitivity, main_geometry->vertex[iMarker], main_geometry->GetnVertex(iMarker), + "MarkerObjectiveCoordinatesTotalSensitivities", true); +} + +vector> CDiscAdjDeformationDriver::GetObjectiveDesignVariablesTotalSensitivities() const { + const auto nDV = GetNumberDesignVariables(); + + vector> values; + + for (auto iDV = 0ul; iDV < nDV; iDV++) { + values.push_back(GetObjectiveDesignVariablesTotalSensitivities(iDV)); + } + + return values; +} + +vector CDiscAdjDeformationDriver::GetObjectiveDesignVariablesTotalSensitivities( + unsigned short iDV) const { + if (iDV >= GetNumberDesignVariables()) { + SU2_MPI::Error("Design Variable index exceeds size.", CURRENT_FUNCTION); + } + + unsigned short nValue = config_container[ZONE_0]->GetnDV_Value(iDV); + vector values(nValue, 0.0); + + for (auto iValue = 0ul; iValue < nValue; iValue++) { + values[iValue] = SU2_TYPE::GetValue(Gradient[iDV][iValue]); + } + + return values; +} diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index cb127ed7d2e..60beb11a56a 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -529,6 +529,19 @@ def main(): test_list.append(pywrapper_wavy_wall_steady) pass_list.append(pywrapper_wavy_wall_steady.run_test()) + # Residual-based discrete adjoint solver + pywrapper_CFD_AD_ResidSolver = TestCase('pywrapper_CFD_AD_ResidSolver') + pywrapper_CFD_AD_ResidSolver.cfg_dir = "py_wrapper/disc_adj_residual_solver" + pywrapper_CFD_AD_ResidSolver.cfg_file = "flow_adjoint.cfg" + pywrapper_CFD_AD_ResidSolver.test_iter = 100 + pywrapper_CFD_AD_ResidSolver.test_vals = [-3.562562, -8.932563, -0.000000, 0.005608] + pywrapper_CFD_AD_ResidSolver.command = TestCase.Command("mpirun -n 2", "python", "run_adjoint.py --parallel -f") + pywrapper_CFD_AD_ResidSolver.timeout = 1600 + pywrapper_CFD_AD_ResidSolver.tol = 0.000001 + pywrapper_CFD_AD_ResidSolver.new_output = False + test_list.append(pywrapper_CFD_AD_ResidSolver) + pass_list.append(pywrapper_CFD_AD_ResidSolver.run_test()) + #################################################################### ### Unsteady Disc. adj. compressible RANS restart optimization ### #################################################################### diff --git a/TestCases/py_wrapper/disc_adj_residual_solver/airfoil/flow_adjoint.cfg b/TestCases/py_wrapper/disc_adj_residual_solver/airfoil/flow_adjoint.cfg new file mode 100644 index 00000000000..dd46f573920 --- /dev/null +++ b/TestCases/py_wrapper/disc_adj_residual_solver/airfoil/flow_adjoint.cfg @@ -0,0 +1,124 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Test residual-based adjoint solver with Python wrapper. % +% Author: H. Patel % +% Date: 11 February 2025 % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER = EULER +KIND_TURB_MODEL = NONE +KIND_DISC_ADJ = RESIDUALS +MATH_PROBLEM = DISCRETE_ADJOINT +OBJECTIVE_FUNCTION = LIFT +ITER = 9999 + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER = 0.3 +AOA = 5.0 +REF_DIMENSIONALIZATION = DIMENSIONAL +FREESTREAM_TEMPERATURE = 293.0 +REYNOLDS_NUMBER = 7E6 +REYNOLDS_LENGTH = 1.0 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH = 1.0 +REF_AREA = 1.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_EULER = ( airfoil ) +MARKER_FAR = ( farfield ) +MARKER_DEFORM_MESH = ( airfoil ) +MARKER_FLUID_LOAD = ( airfoil ) +MARKER_PLOTTING = ( airfoil ) +MARKER_MONITORING = ( airfoil ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD = WEIGHTED_LEAST_SQUARES +CFL_NUMBER = 100000.0 +CFL_ADAPT = NO +CFL_ADAPT_PARAM = ( 0.66667, 2.0, 1.0, 10000.0 ) +RK_ALPHA_COEFF = ( 0.66667, 0.66667, 1.000000 ) + +LINEAR_SOLVER = FGMRES +LINEAR_SOLVER_PREC = ILU +LINEAR_SOLVER_ERROR = 1E-12 +LINEAR_SOLVER_ITER = 1000 +QUASI_NEWTON_NUM_SAMPLES = 10 + +MGLEVEL = 0 +MGCYCLE = V_CYCLE + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW = JST +MUSCL_FLOW = NO +JST_SENSOR_COEFF = ( 0.5, 0.01 ) +TIME_DISCRE_FLOW = EULER_IMPLICIT + +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB = SCALAR_UPWIND +MUSCL_TURB = NO +TIME_DISCRE_TURB = EULER_IMPLICIT + +% --------------------------- MESH PARAMETERS ---------------------------------% +% +DEFORM_MESH = YES +DEFORM_STIFFNESS_TYPE = WALL_DISTANCE +DEFORM_POISSONS_RATIO = 1E6 + +DEFORM_LINEAR_SOLVER = CONJUGATE_GRADIENT +DEFORM_LINEAR_SOLVER_PREC = ILU +DEFORM_LINEAR_SOLVER_ERROR = 1E-8 +DEFORM_LINEAR_SOLVER_ITER = 1000 +DEFORM_CONSOLE_OUTPUT = NO + +READ_BINARY_RESTART = NO + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD = REL_RMS_DENSITY +CONV_RESIDUAL_MINVAL = -12 +WINDOW_CAUCHY_CRIT = NO +CONV_WINDOW_FIELD = (TAVG_DRAG, TAVG_LIFT) +CONV_WINDOW_STARTITER = 0 +CONV_WINDOW_CAUCHY_EPS = 1E-3 +CONV_WINDOW_CAUCHY_ELEMS = 10 + +WINDOW_START_ITER= 1000 +WINDOW_FUNCTION= HANN_SQUARE + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +HISTORY_WRT_FREQ_INNER = 1 +SCREEN_WRT_FREQ_INNER = 1 + +MESH_FILENAME = naca_mesh.su2 +MESH_FORMAT = SU2 + +MESH_OUT_FILENAME = mesh_deformed.su2 +SOLUTION_FILENAME = solution_flow.dat +RESTART_FILENAME = restart_flow.dat +RESTART_ADJ_FILENAME = restart_adjoint_flow.dat +SOLUTION_ADJ_FILENAME = solution_adjoint_flow.dat + +TABULAR_FORMAT = TECPLOT +OUTPUT_FILES = (RESTART, TECPLOT, TECPLOT_ASCII, SURFACE_TECPLOT, SURFACE_TECPLOT_ASCII, CSV, SURFACE_CSV) + +CONV_FILENAME = history +VOLUME_FILENAME = flow +SURFACE_FILENAME = surface_flow + +SCREEN_OUTPUT = (INNER_ITER, RMS_DENSITY, REL_RMS_DENSITY, AVG_CFL, RMS_ADJ_DENSITY) +HISTORY_OUTPUT = ITER, FLOW_COEFF, AERO_COEFF, AOA, RMS_RES, REL_RMS_RES, SENSITIVITY diff --git a/TestCases/py_wrapper/disc_adj_residual_solver/airfoil/flow_forward.cfg b/TestCases/py_wrapper/disc_adj_residual_solver/airfoil/flow_forward.cfg new file mode 100644 index 00000000000..b9568db4c88 --- /dev/null +++ b/TestCases/py_wrapper/disc_adj_residual_solver/airfoil/flow_forward.cfg @@ -0,0 +1,123 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Test residual-based adjoint solver with Python wrapper. % +% Author: H. Patel % +% Date: 11 February 2025 % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER = EULER +KIND_TURB_MODEL = NONE +KIND_DISC_ADJ = RESIDUALS +MATH_PROBLEM = DIRECT +OBJECTIVE_FUNCTION = LIFT +ITER = 9999 + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER = 0.3 +AOA = 5.0 +REF_DIMENSIONALIZATION = DIMENSIONAL +FREESTREAM_TEMPERATURE = 293.0 +REYNOLDS_NUMBER = 7E6 +REYNOLDS_LENGTH = 1.0 + +% ---------------------- REFERENCE VALUE DEFINITION ---------------------------% +% +REF_ORIGIN_MOMENT_X = 0.25 +REF_ORIGIN_MOMENT_Y = 0.00 +REF_ORIGIN_MOMENT_Z = 0.00 +REF_LENGTH = 1.0 +REF_AREA = 1.0 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_EULER = ( airfoil ) +MARKER_FAR = ( farfield ) +MARKER_DEFORM_MESH = ( airfoil ) +MARKER_FLUID_LOAD = ( airfoil ) +MARKER_PLOTTING = ( airfoil ) +MARKER_MONITORING = ( airfoil ) + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD = WEIGHTED_LEAST_SQUARES +CFL_NUMBER = 100.0 +CFL_ADAPT = NO +CFL_ADAPT_PARAM = ( 0.66667, 2.0, 1.0, 10000.0 ) +RK_ALPHA_COEFF = ( 0.66667, 0.66667, 1.000000 ) + +LINEAR_SOLVER = FGMRES +LINEAR_SOLVER_PREC = ILU +LINEAR_SOLVER_ERROR = 1E-5 +LINEAR_SOLVER_ITER = 5 + +MGLEVEL = 0 +MGCYCLE = V_CYCLE + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW = JST +MUSCL_FLOW = NO +JST_SENSOR_COEFF = ( 0.5, 0.01 ) +TIME_DISCRE_FLOW = EULER_IMPLICIT + +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB = SCALAR_UPWIND +MUSCL_TURB = NO +TIME_DISCRE_TURB = EULER_IMPLICIT + +% --------------------------- MESH PARAMETERS ---------------------------------% +% +DEFORM_MESH = YES +DEFORM_STIFFNESS_TYPE = WALL_DISTANCE +DEFORM_POISSONS_RATIO = 1E6 + +DEFORM_LINEAR_SOLVER = CONJUGATE_GRADIENT +DEFORM_LINEAR_SOLVER_PREC = ILU +DEFORM_LINEAR_SOLVER_ERROR = 1E-8 +DEFORM_LINEAR_SOLVER_ITER = 1000 +DEFORM_CONSOLE_OUTPUT = NO + +READ_BINARY_RESTART = NO + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +CONV_FIELD = REL_RMS_DENSITY +CONV_RESIDUAL_MINVAL = -12 +WINDOW_CAUCHY_CRIT = NO +CONV_WINDOW_FIELD = (TAVG_DRAG, TAVG_LIFT) +CONV_WINDOW_STARTITER = 0 +CONV_WINDOW_CAUCHY_EPS = 1E-3 +CONV_WINDOW_CAUCHY_ELEMS = 10 + +WINDOW_START_ITER= 1000 +WINDOW_FUNCTION= HANN_SQUARE + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +HISTORY_WRT_FREQ_INNER = 1 +SCREEN_WRT_FREQ_INNER = 1 + +MESH_FILENAME = naca_mesh.su2 +MESH_FORMAT = SU2 + +MESH_OUT_FILENAME = mesh_deformed.su2 +SOLUTION_FILENAME = solution_flow.dat +RESTART_FILENAME = restart_flow.dat +RESTART_ADJ_FILENAME = restart_adjoint_flow.dat +SOLUTION_ADJ_FILENAME = solution_adjoint_flow.dat + +TABULAR_FORMAT = TECPLOT +OUTPUT_FILES = (RESTART, TECPLOT, TECPLOT_ASCII, SURFACE_TECPLOT, SURFACE_TECPLOT_ASCII, CSV, SURFACE_CSV) + +CONV_FILENAME = history +VOLUME_FILENAME = flow +SURFACE_FILENAME = surface_flow + +SCREEN_OUTPUT = (INNER_ITER, RMS_DENSITY, REL_RMS_DENSITY, AVG_CFL, LIFT, DRAG) +HISTORY_OUTPUT = ITER, FLOW_COEFF, AERO_COEFF, RMS_RES, REL_RMS_RES diff --git a/TestCases/py_wrapper/disc_adj_residual_solver/box/flow_adjoint.cfg b/TestCases/py_wrapper/disc_adj_residual_solver/box/flow_adjoint.cfg new file mode 100644 index 00000000000..b6ef23cec7b --- /dev/null +++ b/TestCases/py_wrapper/disc_adj_residual_solver/box/flow_adjoint.cfg @@ -0,0 +1,111 @@ +SOLVER= EULER +DEFORM_MESH= NO +KIND_TURB_MODEL= SST +MATH_PROBLEM= DISCRETE_ADJOINT +KIND_DISC_ADJ= RESIDUALS +AXISYMMETRIC= NO +RESTART_SOL= NO +SYSTEM_MEASUREMENTS= SI +MGLEVEL= 0 +FLUID_MODEL= IDEAL_GAS +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 287.058 +MACH_NUMBER= 0.25 +REYNOLDS_NUMBER= 1000000 +AOA= 1.0 +SIDESLIP_ANGLE= 0.0 +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_PRESSURE= 101325 +FREESTREAM_TEMPERATURE= 300 +FREESTREAM_DENSITY= 1.24 +REF_ORIGIN_MOMENT_X= 0.00 +REF_ORIGIN_MOMENT_Y= 0.00 +REF_ORIGIN_MOMENT_Z= 0.00 +REF_LENGTH= 1.0 +REF_AREA= 0 +REF_DIMENSIONALIZATION= DIMENSIONAL +MARKER_EULER= ( bottom, top ) +% MARKER_HEATFLUX= ( bottom, top ) +MARKER_FAR= ( left, right ) +MARKER_FLUID_LOAD= ( bottom, top ) +MARKER_DEFORM_MESH= ( bottom, top ) +MARKER_PLOTTING= ( left, bottom, right, top ) +MARKER_MONITORING= ( bottom, top ) +MARKER_DESIGNING= ( bottom, top ) +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +TIME_DISCRE_TURB= EULER_EXPLICIT +CFL_REDUCTION_TURB= 1.0 +NUM_METHOD_GRAD= GREEN_GAUSS +OBJECTIVE_FUNCTION= LIFT +CFL_NUMBER= 100.0 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 0.1, 1.2, 1, 150 ) +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +ITER= 20000 +INNER_ITER= 500 +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-10 +LINEAR_SOLVER_ITER= 500 +DEFORM_LINEAR_SOLVER= CONJUGATE_GRADIENT +DEFORM_LINEAR_SOLVER_PREC= ILU +DEFORM_LINEAR_SOLVER_ERROR= 1E-20 +DEFORM_LINEAR_SOLVER_ITER= 1000 +DISCADJ_LIN_SOLVER= FGMRES +DISCADJ_LIN_PREC= ILU +SMOOTH_GRADIENT= NO +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +VENKAT_LIMITER_COEFF= 0.1 +JST_SENSOR_COEFF= ( 0.5, 0.02 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT +CONV_NUM_METHOD_ADJFLOW= ROE +MUSCL_ADJFLOW= YES +SLOPE_LIMITER_ADJFLOW= NONE +CFL_REDUCTION_ADJFLOW= 1.0 +TIME_DISCRE_ADJFLOW= EULER_EXPLICIT +CONV_FIELD= REL_RMS_DENSITY, LIFT +CONV_RESIDUAL_MINVAL= -16 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-3 +SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, REL_RMS_DENSITY, AVG_CFL, RMS_ADJ_DENSITY) +HISTORY_OUTPUT= ITER, FLOW_COEFF, AERO_COEFF, AOA, RMS_RES, REL_RMS_RES, SENSITIVITY +VOLUME_OUTPUT= (COORDINATES, SOLUTION, PRIMITIVE, RESIDUAL, SENSITIVITY) +MESH_FILENAME= mesh.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= box.deformed.su2 +SOLUTION_FILENAME= solution.dat +SOLUTION_ADJ_FILENAME= restart_adj.dat +TABULAR_FORMAT= TECPLOT +OUTPUT_FILES= (RESTART, TECPLOT, TECPLOT_ASCII, SURFACE_TECPLOT, SURFACE_TECPLOT_ASCII, CSV, SURFACE_CSV) +CONV_FILENAME= history +VOLUME_FILENAME= flow +VOLUME_ADJ_FILENAME= flow_adjoint +GRAD_OBJFUNC_FILENAME= of_grad.dat +SURFACE_FILENAME= surface_flow +SURFACE_ADJ_FILENAME= surface_adjoint +WRT_FORCES_BREAKDOWN= YES +BREAKDOWN_FILENAME= forces_breakdown.dat +OPT_ITERATIONS= 100 +OPT_ACCURACY= 1e-10 +OPT_RELAX_FACTOR= 1.0 +OPT_GRADIENT_FACTOR= 1.0 +OPT_BOUND_UPPER= 10000000000.0 +OPT_BOUND_LOWER= -10000000000.0 +OPT_COMBINE_OBJECTIVE= NO +OPT_CONSTRAINT= NONE +VALUE_OBJFUNC_FILENAME= of_eval.dat +TARGET_CL= 0.0 +MARKER_OUTLET= (NONE) +MULTIPOINT_WEIGHT= (1.0) +MULTIPOINT_MACH_NUMBER= (0.25) +MULTIPOINT_AOA= (1.0) +MULTIPOINT_SIDESLIP_ANGLE= (0.0) +MULTIPOINT_REYNOLDS_NUMBER= (1000000) +MULTIPOINT_TARGET_CL= (0.0) +MULTIPOINT_FREESTREAM_PRESSURE= (101325) +MULTIPOINT_FREESTREAM_TEMPERATURE= (300) +MULTIPOINT_OUTLET_VALUE= (0.0) +MULTIPOINT_MESH_FILENAME= (box.su2) diff --git a/TestCases/py_wrapper/disc_adj_residual_solver/box/flow_forward.cfg b/TestCases/py_wrapper/disc_adj_residual_solver/box/flow_forward.cfg new file mode 100644 index 00000000000..f112b125992 --- /dev/null +++ b/TestCases/py_wrapper/disc_adj_residual_solver/box/flow_forward.cfg @@ -0,0 +1,106 @@ +SOLVER= EULER +DEFORM_MESH= NO +KIND_TURB_MODEL= SST +MATH_PROBLEM= DIRECT +KIND_DISC_ADJ= RESIDUALS +AXISYMMETRIC= NO +SYSTEM_MEASUREMENTS= SI +MGLEVEL= 0 +FLUID_MODEL= IDEAL_GAS +GAMMA_VALUE= 1.4 +GAS_CONSTANT= 287.058 +MACH_NUMBER= 0.25 +REYNOLDS_NUMBER= 1000000 +AOA= 1.0 +SIDESLIP_ANGLE= 0.0 +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_PRESSURE= 101325 +FREESTREAM_TEMPERATURE= 300 +FREESTREAM_DENSITY= 1.24 +REF_ORIGIN_MOMENT_X= 0.00 +REF_ORIGIN_MOMENT_Y= 0.00 +REF_ORIGIN_MOMENT_Z= 0.00 +REF_LENGTH= 1.0 +REF_AREA= 0 +REF_DIMENSIONALIZATION= DIMENSIONAL +MARKER_EULER= ( bottom, top ) +% MARKER_HEATFLUX= ( bottom, top ) +MARKER_FAR= ( left, right ) +MARKER_FLUID_LOAD= ( bottom, top ) +MARKER_DEFORM_MESH= ( bottom, top ) +MARKER_PLOTTING= ( left, bottom, right, top ) +MARKER_MONITORING= ( bottom, top ) +MARKER_DESIGNING= ( bottom, top ) +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +TIME_DISCRE_TURB= EULER_EXPLICIT +CFL_REDUCTION_TURB= 1.0 +NUM_METHOD_GRAD= GREEN_GAUSS +OBJECTIVE_FUNCTION= LIFT +CFL_NUMBER= 10.0 +CFL_ADAPT= NO +CFL_ADAPT_PARAM= ( 0.1, 1.2, 1, 150 ) +RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 ) +ITER= 1000 +LINEAR_SOLVER= BCGSTAB +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ERROR= 1E-10 +LINEAR_SOLVER_ITER= 10 +DEFORM_LINEAR_SOLVER= CONJUGATE_GRADIENT +DEFORM_LINEAR_SOLVER_PREC= ILU +DEFORM_LINEAR_SOLVER_ERROR= 1E-20 +DEFORM_LINEAR_SOLVER_ITER= 1000 +CONV_NUM_METHOD_FLOW= ROE +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE +VENKAT_LIMITER_COEFF= 0.1 +JST_SENSOR_COEFF= ( 0.5, 0.02 ) +TIME_DISCRE_FLOW= EULER_IMPLICIT +CONV_NUM_METHOD_ADJFLOW= ROE +MUSCL_ADJFLOW= YES +SLOPE_LIMITER_ADJFLOW= NONE +CFL_REDUCTION_ADJFLOW= 1.0 +TIME_DISCRE_ADJFLOW= EULER_EXPLICIT +CONV_FIELD= REL_RMS_DENSITY, LIFT +CONV_RESIDUAL_MINVAL= -13 +CONV_STARTITER= 10 +CONV_CAUCHY_ELEMS= 100 +CONV_CAUCHY_EPS= 1E-3 +SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, REL_RMS_DENSITY, AVG_CFL, RMS_ADJ_DENSITY) +HISTORY_OUTPUT= ITER, FLOW_COEFF, AERO_COEFF, AOA, RMS_RES, REL_RMS_RES, SENSITIVITY +VOLUME_OUTPUT= (COORDINATES, SOLUTION, PRIMITIVE, RESIDUAL, SENSITIVITY) +MESH_FILENAME= mesh.su2 +MESH_FORMAT= SU2 +MESH_OUT_FILENAME= box.deformed.su2 +SOLUTION_FILENAME= solution.dat +SOLUTION_ADJ_FILENAME= solution_adj.dat +TABULAR_FORMAT= TECPLOT +OUTPUT_FILES= (RESTART, TECPLOT, TECPLOT_ASCII, SURFACE_TECPLOT, SURFACE_TECPLOT_ASCII, CSV, SURFACE_CSV) +CONV_FILENAME= history +VOLUME_FILENAME= flow +VOLUME_ADJ_FILENAME= flow_adjoint +GRAD_OBJFUNC_FILENAME= of_grad.dat +SURFACE_FILENAME= surface_flow +SURFACE_ADJ_FILENAME= surface_adjoint +WRT_FORCES_BREAKDOWN= YES +BREAKDOWN_FILENAME= forces_breakdown.dat +OPT_ITERATIONS= 100 +OPT_ACCURACY= 1e-10 +OPT_RELAX_FACTOR= 1.0 +OPT_GRADIENT_FACTOR= 1.0 +OPT_BOUND_UPPER= 10000000000.0 +OPT_BOUND_LOWER= -10000000000.0 +OPT_COMBINE_OBJECTIVE= NO +OPT_CONSTRAINT= NONE +VALUE_OBJFUNC_FILENAME= of_eval.dat +TARGET_CL= 0.0 +MARKER_OUTLET= (NONE) +MULTIPOINT_WEIGHT= (1.0) +MULTIPOINT_MACH_NUMBER= (0.25) +MULTIPOINT_AOA= (1.0) +MULTIPOINT_SIDESLIP_ANGLE= (0.0) +MULTIPOINT_REYNOLDS_NUMBER= (1000000) +MULTIPOINT_TARGET_CL= (0.0) +MULTIPOINT_FREESTREAM_PRESSURE= (101325) +MULTIPOINT_FREESTREAM_TEMPERATURE= (300) +MULTIPOINT_OUTLET_VALUE= (0.0) +MULTIPOINT_MESH_FILENAME= (box.su2) diff --git a/TestCases/py_wrapper/disc_adj_residual_solver/run_adjoint.py b/TestCases/py_wrapper/disc_adj_residual_solver/run_adjoint.py new file mode 100644 index 00000000000..8eac3a88db4 --- /dev/null +++ b/TestCases/py_wrapper/disc_adj_residual_solver/run_adjoint.py @@ -0,0 +1,96 @@ +import sys +import pysu2ad +from optparse import OptionParser +import numpy as np + +def matrix_to_numpy(mat): + """Convert a CPyWrapperMatrixView to a numpy array.""" + num_rows, num_cols = mat.Shape() + arr = np.zeros((num_rows, num_cols)) + for i in range(num_rows): + row = mat.Get(i) + for j in range(num_cols): + arr[i, j] = row[j] + + return arr + + +def main(): + # command line options + parser = OptionParser() + parser.add_option("-f", "--file", dest="filename", help="Read config from FILE", metavar="FILE") + parser.add_option("--parallel", action="store_true", help="Specify if we need to initialize MPI", dest="with_MPI", + default=False) + + (options, args) = parser.parse_args() + options.nDim = 2 + options.nZone = 1 + + # initialize MPI + if options.with_MPI: + from mpi4py import MPI + comm = MPI.COMM_WORLD + rank = MPI.COMM_WORLD.Get_rank() + else: + comm = 0 + rank = 0 + + # initialize and preprocess solver + adj_driver = pysu2ad.CDiscAdjSinglezoneDriver(options.filename, options.nZone, comm) + + sys.stdout.flush() + if options.with_MPI: + comm.Barrier() + + adj_driver.Preprocess(0) + adj_driver.Run() + adj_driver.Postprocess() + adj_driver.Update() + adj_driver.Monitor(0) + adj_driver.Output(0) + + # extract sensitivities + if rank == 0: + print("\n------------------------------ Sensitivities -----------------------------\n") + + # farfield variable sensitivities + obj_farfield = adj_driver.GetObjectiveFarfieldVariablesSensitivities() + + if rank == 0: + print("Farfield variable sensitivities (dObjective/dVariables):") + print(f" dI/dMach = {obj_farfield[0]:.6e}") + print(f" dI/dAoA = {obj_farfield[1]:.6e}") + + # dI/dq (nPoint, nVar) + obj_states = adj_driver.ObjectiveStatesSensitivities() + obj_states_np = matrix_to_numpy(obj_states) + if rank == 0: + print(f"\ndObjective/dStates: shape = {obj_states_np.shape}") + print(obj_states_np) + + # dA/dq^T * psi (nPoint, nVar) + res_states = adj_driver.ResidualsStatesSensitivities() + res_states_np = matrix_to_numpy(res_states) + if rank == 0: + print(f"\ndResiduals/dStates^T * psi: shape = {res_states_np.shape}") + print(res_states_np) + + # dI/dxv (nPoint, nDim) + obj_coords = adj_driver.ObjectiveCoordinatesSensitivities() + obj_coords_np = matrix_to_numpy(obj_coords) + if rank == 0: + print(f"\ndObjective/dCoordinates: shape = {obj_coords_np.shape}") + print(obj_coords_np) + + # dA/dxv^T * psi (nPoint, nDim) + res_coords = adj_driver.ResidualsCoordinatesSensitivities() + res_coords_np = matrix_to_numpy(res_coords) + if rank == 0: + print(f"\ndResiduals/dCoordinates^T * psi: shape = {res_coords_np.shape}") + print(res_coords_np) + + adj_driver.Finalize() + + +if __name__ == '__main__': + main() diff --git a/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py b/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py index 286f170b044..af21cbcd44b 100755 --- a/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py +++ b/TestCases/py_wrapper/flatPlate_rigidMotion/launch_flatPlate_rigidMotion.py @@ -76,7 +76,7 @@ def main(): MovingMarker = 'plate' #specified by the user # Get all the tags with the moving option - MovingMarkerList = SU2Driver.GetMarkerTags() + MovingMarkerList = SU2Driver.GetMarkerTags() # Get all the markers defined on this rank and their associated indices. allMarkerIDs = SU2Driver.GetMarkerIndices() diff --git a/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py b/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py index 8ceff23fc30..ca59bcb1381 100755 --- a/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py +++ b/TestCases/py_wrapper/flatPlate_unsteady_CHT/launch_unsteady_CHT_FlatPlate.py @@ -75,7 +75,7 @@ def main(): CHTMarker = 'plate' # Specified by the user # Get all the tags with the CHT option - CHTMarkerList = SU2Driver.GetCHTMarkerTags() + CHTMarkerList = SU2Driver.GetCHTMarkerTags() # Get all the markers defined on this rank and their associated indices. allMarkerIDs = SU2Driver.GetMarkerIndices()