Merge branch 'develop' into develop

Author: pcarruscag

Common/include/CConfig.hpp

@@ -8658,7 +8658,6 @@ class CConfig {
   unsigned short GetKind_Average(void) const { return Kind_Average; }
 
   /*!
-   *
    * \brief Get the direct differentation method.
    * \return direct differentiation method.
    */

SU2_CFD/include/drivers/CDriverBase.hpp

@@ -473,6 +473,25 @@ class CDriverBase {
         "MarkerSolutionTimeN of " + solver->GetSolverName(), false);
   }
 
+  /*!
+   * \brief Get a read/write view of the solution at time N-1 on all mesh nodes of a solver.
+   */
+  inline CPyWrapperMatrixView SolutionTimeN1(unsigned short iSolver) {
+    auto* solver = GetSolverAndCheckMarker(iSolver);
+    return CPyWrapperMatrixView(
+        solver->GetNodes()->GetSolution_time_n1(), "SolutionTimeN1 of " + solver->GetSolverName(), false);
+  }
+
+  /*!
+   * \brief Get a read/write view of the solution at time N-1 on the mesh nodes of a marker.
+   */
+  inline CPyWrapperMarkerMatrixView MarkerSolutionTimeN1(unsigned short iSolver, unsigned short iMarker) {
+    auto* solver = GetSolverAndCheckMarker(iSolver, iMarker);
+    return CPyWrapperMarkerMatrixView(
+        solver->GetNodes()->GetSolution_time_n1(), main_geometry->vertex[iMarker], main_geometry->GetnVertex(iMarker),
+        "MarkerSolutionTimeN1 of " + solver->GetSolverName(), false);
+  }
+
   /*!
    * \brief Get the flow solver primitive variable names with their associated indices.
    * These correspond to the column indices in the matrix returned by Primitives.

SU2_CFD/include/fluid/CCoolProp.hpp

@@ -48,6 +48,7 @@ class CCoolProp final : public CFluidModel {
   su2double Temperature_Critical{0.0}; /*!< \brief critical temperature */
   su2double acentric_factor{0.0};      /*!< \brief acentric factor */
   const su2double dp{0.01};            /*!< threshold for pressure */
+  const su2double dt{0.01};            /*!< threshold for temperature */
 #ifdef USE_COOLPROP
   std::unique_ptr<CoolProp::AbstractState> fluid_entity; /*!< \brief fluid entity */
 #endif
@@ -62,6 +63,18 @@ class CCoolProp final : public CFluidModel {
       Pressure = fmin(Pressure, (1 - dp) * Pressure_Critical);
   }
 
+  /*!
+   * \brief Avoid critical temperature
+   * \param[in,out] Temperature: Modified so that it is not too close to critical temperature to avoid issues in
+   * CoolProp.
+   */
+  void CheckTemperature(su2double& Temperature) const {
+    if (Temperature > Temperature_Critical)
+      Temperature = fmax(Temperature, (1 + dt) * Temperature_Critical);
+    else
+      Temperature = fmin(Temperature, (1 + dt) * Temperature_Critical);
+  }
+
  public:
   /*!
    * \brief Constructor of the class.
@@ -155,4 +168,4 @@ class CCoolProp final : public CFluidModel {
    * \return Value of the constant: Gamma
    */
   su2double GetGamma(void) const { return Gamma; }
-};
+};

SU2_CFD/include/output/CFlowOutput.hpp

@@ -215,7 +215,8 @@ class CFlowOutput : public CFVMOutput{
    * \brief Helper for custom outputs, converts variable names to indices and pointers which are then used
    * to evaluate the custom expressions.
    */
-  void ConvertVariableSymbolsToIndices(const CPrimitiveIndices<unsigned long>& idx, CustomOutput& output) const;
+  void ConvertVariableSymbolsToIndices(const CPrimitiveIndices<unsigned long>& idx, bool allowSkip,
+                                       CustomOutput& output) const;
 
   /*!
    * \brief Compute value of the Q criteration for vortex idenfitication
@@ -296,7 +297,7 @@ class CFlowOutput : public CFVMOutput{
    * \param[in] force_writing - boolean that forces writing of volume output
    * \param[in] iFile - index to the file that we need to consider for volume output
    */
-  bool WriteVolumeOutput(CConfig *config, unsigned long Iter, bool force_writing, unsigned short iFile) override; 
+  bool WriteVolumeOutput(CConfig *config, unsigned long Iter, bool force_writing, unsigned short iFile) override;
   /*!
    * \brief Write the forces breakdown file
    * \param[in] config - Definition of the particular problem per zone.

SU2_CFD/include/output/COutput.hpp

@@ -223,6 +223,10 @@ class COutput {
      We store pointers to the required outputs to speed-up access. ---*/
     std::vector<const su2double*> otherOutputs;
 
+    /*--- For discrete adjoint we may need to skip some expressions because there is one output class
+     for the primal solver and one for the discrete adjoint (each with different variables). ---*/
+    bool skip = false;
+
     /*--- For evaluation, "vars" is a functor (i.e. has operator()) that returns the value of a variable at a given
      point. For example, it can be a wrapper to the primitives pointer, in which case varIndices needs to be setup
      with primitive indices. ---*/
@@ -807,6 +811,14 @@ class COutput {
    */
   void SetCustomOutputs(const CConfig *config);
 
+  /*!
+   * \brief Evaluates function-type custom outputs.
+   * Derived classes can use this to compute simple expressions of other outputs if they
+   * do not implement surface averages. This should be called just before evaluating the
+   * custom objective function.
+   */
+  void ComputeSimpleCustomOutputs(const CConfig *config);
+
   /*!
    * \brief Load values of the history fields common for all solvers.
    * \param[in] config - Definition of the particular problem.

SU2_CFD/include/variables/CVariable.hpp

@@ -500,6 +500,7 @@ class CVariable {
    * \return Pointer to the solution (at time n-1) vector.
    */
   inline su2double *GetSolution_time_n1(unsigned long iPoint) { return Solution_time_n1[iPoint]; }
+  inline MatrixType& GetSolution_time_n1() { return Solution_time_n1; }
 
   /*!
    * \brief Set the value of the old residual.

SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp

@@ -129,6 +129,10 @@ CDiscAdjSinglezoneDriver::~CDiscAdjSinglezoneDriver() {
 
 void CDiscAdjSinglezoneDriver::Preprocess(unsigned long TimeIter) {
 
+  /*--- Set the current time iteration in the config and also in the driver
+   * because the python interface doesn't offer an explicit way of doing it. ---*/
+
+  this->TimeIter = TimeIter;
   config_container[ZONE_0]->SetTimeIter(TimeIter);
 
   /*--- Preprocess the adjoint iteration ---*/

SU2_CFD/src/drivers/CSinglezoneDriver.cpp

@@ -110,8 +110,10 @@ void CSinglezoneDriver::StartSolver() {
 
 void CSinglezoneDriver::Preprocess(unsigned long TimeIter) {
 
-  /*--- Set the current time iteration in the config ---*/
+  /*--- Set the current time iteration in the config and also in the driver
+   * because the python interface doesn't offer an explicit way of doing it. ---*/
 
+  this->TimeIter = TimeIter;
   config_container[ZONE_0]->SetTimeIter(TimeIter);
 
   /*--- Store the current physical time in the config container, as

SU2_CFD/src/fluid/CCoolProp.cpp

@@ -56,25 +56,20 @@ void CCoolProp::SetTDState_rhoe(su2double rho, su2double e) {
   dTdrho_e = fluid_entity->first_partial_deriv(CoolProp::iT, CoolProp::iDmass, CoolProp::iUmass);
   dTde_rho = fluid_entity->first_partial_deriv(CoolProp::iT, CoolProp::iUmass, CoolProp::iDmass);
   if (fluid_entity->phase() == CoolProp::iphase_twophase) {
-    // assume it is pure gas
-    fluid_entity->specify_phase(CoolProp::iphase_gas);
+    // impose gas phase
+    Temperature = Temperature + 0.1;
     CheckPressure(Pressure);
+    CheckTemperature(Temperature);
     fluid_entity->update(CoolProp::PT_INPUTS, Pressure, Temperature);
-    if (abs(fluid_entity->rhomass() / Density - 1) < dp) {
-      // origial phase is near saturation gas, then just compute sound speed
-      SoundSpeed2 = pow(fluid_entity->speed_sound(), 2);
-    } else {
-      // original phase is not near saturation gas, then specify the phase as gas phase
-      fluid_entity->specify_phase(CoolProp::iphase_gas);
-      SetTDState_PT(Pressure, Temperature);
-    }
+    SoundSpeed2 = pow(fluid_entity->speed_sound(), 2);
   } else {
     SoundSpeed2 = pow(fluid_entity->speed_sound(), 2);
   }
 }
 
 void CCoolProp::SetTDState_PT(su2double P, su2double T) {
   CheckPressure(P);
+  CheckTemperature(T);
   fluid_entity->update(CoolProp::PT_INPUTS, P, T);
   su2double rho = fluid_entity->rhomass();
   su2double e = fluid_entity->umass();
@@ -130,4 +125,4 @@ CCoolProp::CCoolProp(const string& fluidname) {
       "or autodiff",
       CURRENT_FUNCTION);
 }
-#endif
+#endif