Fix adjoint for streamwise periodic massflow + General handling of adjoints of additional solution variables

Common/include/option_structure.hpp

@@ -2254,6 +2254,8 @@ struct StreamwisePeriodicValues {
   su2double Streamwise_Periodic_MassFlow;           /*!< \brief Value of current massflow [kg/s] which results in a delta p and therefore an artificial body force vector. */
   su2double Streamwise_Periodic_IntegratedHeatFlow; /*!< \brief Value of of the net sum of heatflow [W] into the domain. */
   su2double Streamwise_Periodic_InletTemperature;   /*!< \brief Area avg static Temp [K] at the periodic inlet. Used for adaptive outlet heatsink. */
+  su2double Streamwise_Periodic_BoundaryArea;       /*!< \brief Global Surface area of the streamwise periodic interface. */
+  su2double Streamwise_Periodic_AvgDensity;         /*!< \brief Area avg density on the periodic interface. */
 };
 
 /*!

Common/src/CConfig.cpp

@@ -4803,8 +4803,6 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("A MARKER_PERIODIC pair has to be set with KIND_STREAMWISE_PERIODIC != NONE.", CURRENT_FUNCTION);
     if (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal != 0)
       SU2_MPI::Error("No MARKER_ISOTHERMAL marker allowed with STREAMWISE_PERIODIC_TEMPERATURE= YES, only MARKER_HEATFLUX & MARKER_SYM.", CURRENT_FUNCTION);
-    if (DiscreteAdjoint && Kind_Streamwise_Periodic == ENUM_STREAMWISE_PERIODIC::MASSFLOW)
-      SU2_MPI::Error("Discrete Adjoint currently not validated for prescribed MASSFLOW.", CURRENT_FUNCTION);
     if (Ref_Inc_NonDim != DIMENSIONAL)
       SU2_MPI::Error("Streamwise Periodicity only works with \"INC_NONDIM= DIMENSIONAL\", the nondimensionalization with source terms doesn;t work in general.", CURRENT_FUNCTION);
     if (Axisymmetric)

SU2_CFD/include/solvers/CIncEulerSolver.hpp

@@ -39,7 +39,7 @@
 class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME::INCOMPRESSIBLE> {
 protected:
   vector<CFluidModel*> FluidModel;   /*!< \brief fluid model used in the solver. */
-  StreamwisePeriodicValues SPvals;
+  StreamwisePeriodicValues SPvals, SPvalsUpdated;
 
   /*!
    * \brief Preprocessing actions common to the Euler and NS solvers.
@@ -395,8 +395,30 @@ class CIncEulerSolver : public CFVMFlowSolverBase<CIncEulerVariable, ENUM_REGIME
   inline bool GetHasHybridParallel() const final { return true; }
 
   /*!
-   * \brief Get values for streamwise periodic flow: delta P, m_dot, inlet T, integrated heat.
-   * \return Struct holding 4 su2doubles.
+   * \brief Get values for streamwise periodic flow: delta P, m_dot, inlet T, integrated heat, etc.
+   * \return Struct holding streamwise periodic values.
    */
   StreamwisePeriodicValues GetStreamwisePeriodicValues() const final { return SPvals; }
+
+  /*!
+   * \brief Register In- or Output.
+   * \param[in] input - Boolean whether In- or Output should be registered.
+   * \param[in] config - The particular config.
+   * \returns The number of extra variables.
+   */
+  unsigned long RegisterSolutionExtra(bool input, const CConfig* config) final;
+
+  /*!
+   * \brief Seed the adjoint of the extra solution at the output.
+   * \param[in] adj_sol - Vector containing the adjoint solution to seed.
+   * \param[in] config - The particular config.
+   */
+  void SetAdjoint_SolutionExtra(const su2activevector& adj_sol, const CConfig* config) final;
+
+  /*!
+   * \brief Extract the adjoint of the extra solution at the input.
+   * \param[out] adj_sol - Vector to store the adjoint into.
+   * \param[in] config - The particular config.
+   */
+  void ExtractAdjoint_SolutionExtra(su2activevector& adj_sol, const CConfig* config) final;
 };

SU2_CFD/include/solvers/CSolver.hpp

@@ -3543,6 +3543,28 @@ class CSolver {
    */
   inline virtual void ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm){}
 
+  /*!
+   * \brief Register In- or Output.
+   * \param[in] input - Boolean whether In- or Output should be registered.
+   * \param[in] config - The particular config.
+   * \returns The number of extra variables.
+   */
+  virtual unsigned long RegisterSolutionExtra(bool input, const CConfig* config) { return 0; }
+
+  /*!
+   * \brief Seed the adjoint of the extra solution at the output.
+   * \param[in] adj_sol - Vector containing the adjoint solution to seed.
+   * \param[in] config - The particular config.
+   */
+  virtual void SetAdjoint_SolutionExtra(const su2activevector& adj_sol, const CConfig* config) {}
+
+  /*!
+   * \brief Extract the adjoint of the extra solution at the input.
+   * \param[out] adj_sol - Vector to store the adjoint into.
+   * \param[in] config - The particular config.
+   */
+  virtual void ExtractAdjoint_SolutionExtra(su2activevector& adj_sol, const CConfig* config) {}
+
   /*!
    * \brief  A virtual member.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -4281,8 +4303,8 @@ class CSolver {
   inline virtual bool GetHasHybridParallel() const { return false; }
 
   /*!
-   * \brief Get values for streamwise periodc flow: delta P, m_dot, inlet T, integrated heat.
-   * \return Struct holding 4 su2doubles.
+   * \brief Get values for streamwise periodic flow: delta P, m_dot, inlet T, integrated heat, etc.
+   * \return Struct holding streamwise periodic values.
    */
   virtual StreamwisePeriodicValues GetStreamwisePeriodicValues() const { return StreamwisePeriodicValues(); }
 

SU2_CFD/include/variables/CDiscAdjVariable.hpp

@@ -54,9 +54,10 @@ class CDiscAdjVariable : public CVariable {
   CDiscAdjVariable(const su2double* sol, unsigned long npoint, unsigned long ndim, unsigned long nvar, CConfig *config);
 
   /*!
-   * \brief Destructor of the class.
+   * \brief Allocate extra adjoint variables.
+   * \param[in] nVarExtra - Number of extra variables.
    */
-  ~CDiscAdjVariable() override = default;
+  void AllocateAdjointSolutionExtra(unsigned long nVarExtra);
 
   /*!
    * \brief Set the sensitivity at the node

SU2_CFD/include/variables/CHeatVariable.hpp

@@ -37,8 +37,6 @@
  */
 class CHeatVariable final : public CVariable {
 protected:
-  MatrixType Solution_Direct;  /*!< \brief Direct solution container for use in the adjoint Heat solver. */
-
   CVectorOfMatrix& Gradient_Reconstruction;  /*!< \brief Reference to the gradient of the primitive variables for MUSCL reconstruction for the convective term */
   CVectorOfMatrix Gradient_Aux;              /*!< \brief Auxiliary structure to store a second gradient for reconstruction, if required. */
 

SU2_CFD/include/variables/CVariable.hpp

@@ -94,6 +94,13 @@ class CVariable {
   su2matrix<int> AD_InputIndex;    /*!< \brief Indices of Solution variables in the adjoint vector. */
   su2matrix<int> AD_OutputIndex;   /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */
 
+  VectorType SolutionExtra; /*!< \brief Stores adjoint solution for extra solution variables.
+                                        Currently only streamwise periodic pressure-drop for massflow prescribed flows. */
+  VectorType ExternalExtra; /*!< \brief External storage for the adjoint value (i.e. for the OF mainly */
+
+  VectorType SolutionExtra_BGS_k; /*!< \brief Intermediate storage, enables cross term extraction as that is also pushed to Solution. */
+
+ protected:
   unsigned long nPoint = 0;  /*!< \brief Number of points in the domain. */
   unsigned long nDim = 0;      /*!< \brief Number of dimension of the problem. */
   unsigned long nVar = 0;        /*!< \brief Number of variables of the problem. */
@@ -406,6 +413,31 @@ class CVariable {
     for(unsigned long iVar = 0; iVar < nVar; iVar++) External(iPoint,iVar) += val_sol[iVar];
   }
 
+  /*!
+   * \brief Store the adjoint solution of the extra adjoint into the external container.
+   */
+  void Add_SolutionExtra_To_ExternalExtra() {
+    assert(SolutionExtra.size() == ExternalExtra.size());
+    for (auto iEntry = 0ul; iEntry < SolutionExtra.size(); iEntry++)
+      ExternalExtra[iEntry] = SolutionExtra[iEntry];
+  }
+
+  /*!
+   * \brief Add the external contribution to the solution for the extra adjoint solutions.
+   */
+  void Add_ExternalExtra_To_SolutionExtra() {
+    assert(SolutionExtra.size() == ExternalExtra.size());
+    for (auto iEntry = 0ul; iEntry < SolutionExtra.size(); iEntry++)
+      SolutionExtra[iEntry] += ExternalExtra[iEntry];
+  }
+
+  /*!
+   * \brief Return the extra adjoint solution.
+   * \return Reference to extra adjoint solution.
+   */
+  inline VectorType& GetSolutionExtra() { return SolutionExtra; }
+  inline const VectorType& GetSolutionExtra() const { return SolutionExtra; }
+
   /*!
    * \brief Update the variables using a conservative format.
    * \param[in] iPoint - Point index.

SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp

@@ -430,7 +430,7 @@ void CDiscAdjMultizoneDriver::Run() {
       /*--- Compute residual from Solution and Solution_BGS_k and update the latter. ---*/
 
       SetResidual_BGS(iZone);
-
+      Set_BGSSolution_k_To_Solution(iZone);
     }
 
     /*--- Set the multizone output. ---*/

SU2_CFD/src/drivers/CMultizoneDriver.cpp

@@ -403,8 +403,10 @@ bool CMultizoneDriver::OuterConvergence(unsigned long OuterIter) {
     auto solvers = solver_container[iZone][INST_0][MESH_0];
 
     for (unsigned short iSol = 0; iSol < MAX_SOLS; iSol++){
-      if (solvers[iSol] != nullptr)
+      if (solvers[iSol] != nullptr) {
         solvers[iSol]->ComputeResidual_Multizone(geometry_container[iZone][INST_0][MESH_0], config_container[iZone]);
+        solvers[iSol]->GetNodes()->Set_BGSSolution_k();
+      }
     }
 
     /*--- Make sure that everything is loaded into the output container. ---*/

SU2_CFD/src/output/CAdjFlowIncOutput.cpp

@@ -113,6 +113,10 @@ void CAdjFlowIncOutput::SetHistoryOutputFields(CConfig *config){
   /// DESCRIPTION: Maximum residual of the temperature.
   AddHistoryOutput("RMS_ADJ_TEMPERATURE", "rms[A_T]", ScreenOutputFormat::FIXED, "RMS_RES", "Root-mean square residual of the adjoint temperature.", HistoryFieldType::RESIDUAL);
 
+  if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) {
+    AddHistoryOutput("ADJOINT_SOLEXTRA", "Adjoint_SolExtra", ScreenOutputFormat::FIXED, "ADJOINT_SOLEXTRA", "Adjoint value of the first extra Solution.", HistoryFieldType::COEFFICIENT);
+  }
+
   AddHistoryOutputFields_AdjScalarRMS_RES(config);
 
   if (config->AddRadiation()){
@@ -205,6 +209,11 @@ void CAdjFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CS
   if (config->AddRadiation()){
     SetHistoryOutputValue("RMS_ADJ_RAD_ENERGY", log10(adjrad_solver->GetRes_RMS(0)));
   }
+
+  if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) {
+    SetHistoryOutputValue("ADJOINT_SOLEXTRA", adjflow_solver->GetNodes()->GetSolutionExtra()[0]);
+  }
+
   SetHistoryOutputValue("MAX_ADJ_PRESSURE", log10(adjflow_solver->GetRes_Max(0)));
   SetHistoryOutputValue("MAX_ADJ_VELOCITY-X", log10(adjflow_solver->GetRes_Max(1)));
   SetHistoryOutputValue("MAX_ADJ_VELOCITY-Y", log10(adjflow_solver->GetRes_Max(2)));

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

@@ -29,20 +29,19 @@
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
 #include "../../../Common/include/parallelization/omp_structure.hpp"
 
-CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *direct_solver,
+CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *direct_sol,
                                unsigned short Kind_Solver, unsigned short iMesh)  : CSolver() {
 
+  /*-- Store some information about direct solver ---*/
+
+  KindDirect_Solver = Kind_Solver;
+  direct_solver = direct_sol;
+
   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();
@@ -86,6 +85,11 @@ CDiscAdjSolver::CDiscAdjSolver(CGeometry *geometry, CConfig *config, CSolver *di
   nodes = new CDiscAdjVariable(Solution.data(), nPoint, nDim, nVar, config);
   SetBaseClassPointerToNodes();
 
+  /*--- Allocate extra solution variables, if any are in use. ---*/
+
+  const auto nVarExtra = direct_solver->RegisterSolutionExtra(true, config);
+  nodes->AllocateAdjointSolutionExtra(nVarExtra);
+
   switch(KindDirect_Solver){
   case RUNTIME_FLOW_SYS:
     SolverName = "ADJ.FLOW";
@@ -163,6 +167,8 @@ void CDiscAdjSolver::RegisterSolution(CGeometry *geometry, CConfig *config) {
   /*--- Boolean true indicates that an input is registered ---*/
   direct_solver->GetNodes()->RegisterSolution(true);
 
+  direct_solver->RegisterSolutionExtra(true, config);
+
   if (time_n_needed)
     direct_solver->GetNodes()->RegisterSolution_time_n();
 
@@ -298,6 +304,8 @@ void CDiscAdjSolver::RegisterOutput(CGeometry *geometry, CConfig *config) {
   /*--- Register variables as output of the solver iteration. Boolean false indicates that an output is registered ---*/
 
   direct_solver->GetNodes()->RegisterSolution(false);
+
+  direct_solver->RegisterSolutionExtra(false, config);
 }
 
 void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm) {
@@ -337,6 +345,8 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
   }
   END_SU2_OMP_FOR
 
+  direct_solver->ExtractAdjoint_SolutionExtra(nodes->GetSolutionExtra(), config);
+
   /*--- Residuals and time_n terms are not needed when evaluating multizone cross terms. ---*/
   if (CrossTerm) return;
 
@@ -468,6 +478,8 @@ void CDiscAdjSolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config) {
     direct_solver->GetNodes()->SetAdjointSolution(iPoint,Solution);
   }
   END_SU2_OMP_FOR
+
+  direct_solver->SetAdjoint_SolutionExtra(nodes->GetSolutionExtra(), config);
 }
 
 void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*) {

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -100,6 +100,7 @@ CIncEulerSolver::CIncEulerSolver(CGeometry *geometry, CConfig *config, unsigned
     string filename_ = "flow";
     filename_ = config->GetFilename(filename_, ".meta", Unst_RestartIter);
     Read_SU2_Restart_Metadata(geometry, config, adjoint, filename_);
+    if (rank==MASTER_NODE) cout << "Setting streamwise periodic pressure drop from restart metadata file." << endl;
   }
 
   /*--- Set the gamma value ---*/
@@ -1614,8 +1615,32 @@ void CIncEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_cont
 
     } // if turbulent
 
-    /*--- Set delta_p, m_dot, inlet_T, integrated_heat ---*/
-    numerics->SetStreamwisePeriodicValues(SPvals);
+    if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) {
+      /*---------------------------------------------------------------------------------------------*/
+      /*--- Update the Pressure Drop [Pa] for the Momentum source term if Massflow is prescribed. ---*/
+      /*--- The Pressure drop is iteratively adapted to result in the prescribed Target-Massflow. ---*/
+      /*---------------------------------------------------------------------------------------------*/
+
+      /*--- Compute update to Delta p based on massflow-difference ---*/
+      const su2double Average_Density_Global = SPvals.Streamwise_Periodic_AvgDensity;
+      const su2double Area_Global = SPvals.Streamwise_Periodic_BoundaryArea;
+      const su2double TargetMassFlow = config->GetStreamwise_Periodic_TargetMassFlow() / (config->GetDensity_Ref() * config->GetVelocity_Ref());
+      const su2double MassFlow_Global = SPvals.Streamwise_Periodic_MassFlow;
+      const su2double ddP = 0.5 / ( Average_Density_Global * pow(Area_Global, 2)) * (pow(TargetMassFlow, 2) - pow(MassFlow_Global, 2));
+
+      /*--- Store updated pressure difference ---*/
+      const su2double damping_factor = config->GetInc_Outlet_Damping();
+      SPvalsUpdated = SPvals;
+      SPvalsUpdated.Streamwise_Periodic_PressureDrop += damping_factor*ddP;
+      if (!config->GetDiscrete_Adjoint())
+        SPvals = SPvalsUpdated;
+
+      /*--- Set delta_p, m_dot, inlet_T, integrated_heat ---*/
+      numerics->SetStreamwisePeriodicValues(SPvalsUpdated);
+    }
+    else {
+      numerics->SetStreamwisePeriodicValues(SPvals);
+    }
 
     AD::StartNoSharedReading();
 
@@ -3069,18 +3094,35 @@ void CIncEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf
 }
 
 void CIncEulerSolver::SetFreeStream_Solution(const CConfig *config){
-
-  unsigned long iPoint;
-  unsigned short iDim;
-
   SU2_OMP_FOR_STAT(omp_chunk_size)
-  for (iPoint = 0; iPoint < nPoint; iPoint++){
+  for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){
     nodes->SetSolution(iPoint,0, Pressure_Inf);
-    for (iDim = 0; iDim < nDim; iDim++){
+    for (unsigned short iDim = 0; iDim < nDim; iDim++){
       nodes->SetSolution(iPoint,iDim+1, Velocity_Inf[iDim]);
     }
     nodes->SetSolution(iPoint,nDim+1, Temperature_Inf);
   }
   END_SU2_OMP_FOR
 
 }
+
+unsigned long CIncEulerSolver::RegisterSolutionExtra(bool input, const CConfig* config) {
+  if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) {
+    if (input) AD::RegisterInput(SPvals.Streamwise_Periodic_PressureDrop);
+    else AD::RegisterOutput(SPvalsUpdated.Streamwise_Periodic_PressureDrop);
+    return 1;
+  }
+  return 0;
+}
+
+void CIncEulerSolver::SetAdjoint_SolutionExtra(const su2activevector& adj_sol, const CConfig* config) {
+  if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) {
+    SU2_TYPE::SetDerivative(SPvalsUpdated.Streamwise_Periodic_PressureDrop, SU2_TYPE::GetValue(adj_sol[0]));
+  }
+}
+
+void CIncEulerSolver::ExtractAdjoint_SolutionExtra(su2activevector& adj_sol, const CConfig* config) {
+  if (config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) {
+    adj_sol[0] = SU2_TYPE::GetDerivative(SPvals.Streamwise_Periodic_PressureDrop);
+  }
+}