[WIP] Isothermal streamwise periodicity

Common/include/CConfig.hpp

@@ -1069,7 +1069,9 @@ class CConfig {
   bool Streamwise_Periodic_Temperature;              /*!< \brief Use real periodicity for Energy equation or otherwise outlet source term. */
   su2double Streamwise_Periodic_PressureDrop;        /*!< \brief Value of prescribed pressure drop [Pa] which results in an artificial body force vector. */
   su2double Streamwise_Periodic_TargetMassFlow;      /*!< \brief Value of prescribed massflow [kg/s] which results in an delta p and therefore an artificial body force vector. */
+  su2double Streamwise_Periodic_ComputedMassFlow;      /*!< \brief Value of computed massflow [kg/s] corresponding to a prescribed delta p. */
   su2double Streamwise_Periodic_OutletHeat;          /*!< /brief Heatflux boundary [W/m^2] imposed at streamwise periodic outlet. */
+  su2double Streamwise_Periodic_LambdaL;            /*!< /brief Exp coefficient for iso-thermal BCs Streamwise Periodic. */
 
   su2double *FreeStreamTurboNormal;     /*!< \brief Direction to initialize the flow in turbomachinery computation */
   su2double Restart_Bandwidth_Agg;      /*!< \brief The aggregate of the bandwidth for writing binary restarts (to be averaged later). */
@@ -3187,6 +3189,12 @@ class CConfig {
    */
   unsigned short GetnMarker_HeatFlux(void) const { return nMarker_HeatFlux; }
 
+  /*!
+   * \brief Get the total (local) number of isothermal markers.
+   * \return Total number of isothermal markers.
+   */
+  unsigned short GetnMarker_Isothermal(void) const { return nMarker_Isothermal; }
+
   /*!
    * \brief Get the total number of rough markers.
    * \return Total number of heat flux markers.
@@ -6219,6 +6227,17 @@ class CConfig {
    */
   void SetStreamwise_Periodic_PressureDrop(su2double Streamwise_Periodic_PressureDrop_) { Streamwise_Periodic_PressureDrop = Streamwise_Periodic_PressureDrop_; }
 
+  /*!
+   * \brief Get the value of the MassFlow from which body force vector is computed.
+   * \return MassFlow for body force computation.
+   */
+  su2double GetStreamwise_Periodic_ComputedMassFlow(void) const { return Streamwise_Periodic_ComputedMassFlow; }
+
+  /*!
+   * \brief Set the value of the MassFlow from which body force vector is computed. Necessary for Restart metadata??
+   */
+  void SetStreamwise_Periodic_ComputedMassFlow(su2double Streamwise_Periodic_MassFlow_) { Streamwise_Periodic_ComputedMassFlow = Streamwise_Periodic_MassFlow_; }
+
   /*!
    * \brief Get the value of the massflow from which body force vector is computed.
    * \return Massflow for body force computation.
@@ -9873,4 +9892,16 @@ class CConfig {
    */
   LM_ParsedOptions GetLMParsedOptions() const { return lmParsedOptions; }
 
+  /*!
+   * \brief Set Lambda L for Streamwise Periodic
+   * \return bool
+   */
+  void SetStreamwise_Periodic_LamdaL(su2double value) { Streamwise_Periodic_LambdaL = value; }
+
+  /*!
+   * \brief Get Lambda L for Streamwise Periodic
+   * \return su2double
+   */
+  su2double GetStreamwise_Periodic_LamdaL(void) const { return Streamwise_Periodic_LambdaL; }
+
 };

Common/include/option_structure.hpp

@@ -1993,6 +1993,9 @@ enum ENUM_OBJECTIVE {
   TOPOL_DISCRETENESS = 63,      /*!< \brief Measure of the discreteness of the current topology. */
   TOPOL_COMPLIANCE = 64,        /*!< \brief Measure of the discreteness of the current topology. */
   STRESS_PENALTY = 65,          /*!< \brief Penalty function of VM stresses above a maximum value. */
+  STREAMWISE_LAMBDAL= 71, /*!< /brief temp. expo. coefficient for iso-thermal BCs Streamwise Periodic. */
+  STREAMWISE_PERIODIC_DP= 72, /*!< /brief pressure drop in Streamwise Periodic flow domain. */
+  STREAMWISE_PERIODIC_MASSFLOW= 73, /*!< /brief massflow rate in Streamwise Periodic flow domain. */
 };
 static const MapType<std::string, ENUM_OBJECTIVE> Objective_Map = {
   MakePair("DRAG", DRAG_COEFFICIENT)
@@ -2028,6 +2031,9 @@ static const MapType<std::string, ENUM_OBJECTIVE> Objective_Map = {
   MakePair("SURFACE_PRESSURE_DROP", SURFACE_PRESSURE_DROP)
   MakePair("SURFACE_SPECIES_0", SURFACE_SPECIES_0)
   MakePair("SURFACE_SPECIES_VARIANCE", SURFACE_SPECIES_VARIANCE)
+  MakePair("STREAMWISE_LAMBDAL", STREAMWISE_LAMBDAL)
+  MakePair("STREAMWISE_PERIODIC_DP", STREAMWISE_PERIODIC_DP)
+  MakePair("STREAMWISE_PERIODIC_MASSFLOW", STREAMWISE_PERIODIC_MASSFLOW)
   MakePair("CUSTOM_OBJFUNC", CUSTOM_OBJFUNC)
   MakePair("REFERENCE_GEOMETRY", REFERENCE_GEOMETRY)
   MakePair("REFERENCE_NODE", REFERENCE_NODE)
@@ -2655,6 +2661,8 @@ struct StreamwisePeriodicValues {
   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. */
+  su2double Streamwise_Periodic_LambdaL;            /*!< \brief Temperature Gradient in iso-thermal BCs. */
+  su2double Streamwise_Periodic_ThetaScaling;       /*!< \brief Scaled Temperature for iso-thermal BCs. */
 };
 
 /*!

Common/src/CConfig.cpp

@@ -5150,8 +5150,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       SU2_MPI::Error("Streamwise Periodic Flow + Incompressible Euler: Not tested yet.", CURRENT_FUNCTION);
     if (nMarker_PerBound == 0)
       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 (Energy_Equation && Streamwise_Periodic_Temperature && nMarker_Isothermal > 0 && nMarker_HeatFlux > 0)
+    //   SU2_MPI::Error("MARKER_ISOTHERMAL and MARKER_HEATFLUX are not allowed simultaneously with STREAMWISE_PERIODIC_TEMPERATURE= YES, only one of these is allowed MARKER_ISOTHERMAL, MARKER_HEATFLUX in a fluid zone.", 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/interfaces/CInterface.hpp

@@ -219,4 +219,7 @@ class CInterface {
    */
   void GatherAverageValues(CSolver *donor_solution, CSolver *target_solution, unsigned short donorZone);
 
+  inline virtual void SetLambdaL(CSolver *target_solution, CGeometry *target_geometry, const CConfig *target_config,
+                          unsigned long Marker_Target) {};
+
 };

SU2_CFD/include/interfaces/cht/CConjugateHeatInterface.hpp

@@ -70,4 +70,7 @@ class CConjugateHeatInterface : public CInterface {
    */
   void SetTarget_Variable(CSolver *target_solution, CGeometry *target_geometry, const CConfig *target_config,
                           unsigned long Marker_Target, unsigned long Vertex_Target, unsigned long Point_Target) override;
+
+  void SetLambdaL(CSolver *target_solution, CGeometry *target_geometry, const CConfig *target_config,
+                          unsigned long Marker_Target) override;
 };

SU2_CFD/include/numerics/flow/flow_sources.hpp

@@ -391,6 +391,8 @@ class CSourceIncStreamwise_Periodic final : public CSourceBase_Flow {
   bool turbulent; /*!< \brief Turbulence model used. */
   bool energy;    /*!< \brief Energy equation on. */
   bool streamwisePeriodic_temperature; /*!< \brief Periodicity in energy equation */
+  bool bool_heat_flux_bc; /*!< \brief HEAT Flux BC boolean. */
+  bool bool_isotherml_bc; /*!< \brief ISO THERMAL BC boolean. */
   su2double Streamwise_Coord_Vector[MAXNDIM] = {0.0}; /*!< \brief Translation vector between streamwise periodic surfaces. */
 
   su2double norm2_translation, /*!< \brief Square of distance between the 2 periodic surfaces. */

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -2933,6 +2933,15 @@ su2double CFVMFlowSolverBase<V,R>::EvaluateCommonObjFunc(const CConfig& config)
     case SURFACE_SPECIES_VARIANCE:
       objFun += weight * config.GetSurface_Species_Variance(0);
       break;
+    case STREAMWISE_LAMBDAL:
+      objFun += weight * config.GetStreamwise_Periodic_LamdaL();
+      break;
+    case STREAMWISE_PERIODIC_DP:
+      objFun += weight * config.GetStreamwise_Periodic_PressureDrop();
+      break;
+    case STREAMWISE_PERIODIC_MASSFLOW:
+      objFun += weight * config.GetStreamwise_Periodic_ComputedMassFlow();
+      break;
     case CUSTOM_OBJFUNC:
       objFun += weight * Total_Custom_ObjFunc;
       break;

SU2_CFD/include/solvers/CHeatSolver.hpp

@@ -318,6 +318,20 @@ class CHeatSolver final : public CScalarSolver<CHeatVariable> {
                    CSolver **solver_container,
                    CConfig *config) override;
 
+/*!
+   * \brief Compute the viscous residuals for the turbulent equation.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] solver_container - Container vector with all the solutions.
+   * \param[in] numerics_container - Description of the numerical method.
+   * \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.
+   */
+  void Source_Residual(CGeometry *geometry,
+                        CSolver **solver_container,
+                        CNumerics **numerics_container,
+                        CConfig *config,
+                        unsigned short iMesh) override;
   /*!
    * \brief Get value of the heat load (integrated heat flux).
    * \return Value of the heat load (integrated heat flux).

SU2_CFD/src/interfaces/cht/CConjugateHeatInterface.cpp

@@ -171,3 +171,47 @@
                                               target_config->GetRelaxation_Factor_CHT(), Target_Variable[3]);
   }
 }
+
+void CConjugateHeatInterface::SetLambdaL(CSolver *target_solution, CGeometry *target_geometry,
+                                                 const CConfig *target_config, unsigned long Marker_Target) {
+  unsigned short iPoint, iMarker, iVertex;
+  su2double term1=0.0, term2=0.0, term3=0.0;
+  auto nDim = target_geometry->GetnDim();
+
+    /*--- Loop over all heatflux Markers ---*/
+  for (iMarker = 0; iMarker < target_config->GetnMarker_All(); iMarker++) {
+
+    if (target_config->GetMarker_All_KindBC(iMarker) == CHT_WALL_INTERFACE) {
+
+    /*--- Identify the boundary by string name and retrive heatflux from config ---*/
+    const auto Marker_StringTag = target_config->GetMarker_All_TagBound(iMarker);
+    const su2double Wall_HeatFlux = target_config->GetWall_HeatFlux(Marker_StringTag);
+
+      for (iVertex = 0ul; iVertex < target_geometry->nVertex[iMarker]; iVertex++) {
+
+      auto iPoint = target_geometry->vertex[iMarker][iVertex]->GetNode();
+
+      if (!target_geometry->nodes->GetDomain(iPoint)) continue;
+
+      const auto AreaNormal = target_geometry->vertex[iMarker][iVertex]->GetNormal();
+
+      const su2double FaceArea = GeometryToolbox::Norm(nDim, AreaNormal);
+
+      // compute the constant
+      const auto GradT = target_solution->GetNodes()->GetGradient_Primitive(iPoint)[0]; //TODO: check index
+
+      term1 += GeometryToolbox::DotProduct(nDim, GradT, AreaNormal);
+
+      // compute coefficient of lambda
+      term2 += target_solution->GetNodes()->GetTemperature(iPoint) * FaceArea;
+      }
+    } // loop Vertices
+  }
+
+  // compute coefficient of lambda^2
+for (iPoint = 0; iPoint<=target_geometry->GetnPointDomain(); iPoint++) {
+  term3 += target_solution->GetNodes()->GetTemperature(iPoint) * target_geometry->nodes->GetVolume(iPoint);
+}
+
+  // Compute LambdaL using the quadratic equation
+}

SU2_CFD/src/numerics/flow/flow_sources.cpp

@@ -723,6 +723,8 @@
   turbulent = (config->GetKind_Turb_Model() != TURB_MODEL::NONE);
   energy    = config->GetEnergy_Equation();
   streamwisePeriodic_temperature = config->GetStreamwise_Periodic_Temperature();
+  bool_heat_flux_bc = (config->GetnMarker_HeatFlux() > 0);
+  bool_isotherml_bc = (config->GetnMarker_Isothermal() > 0);
 
   for (unsigned short iDim = 0; iDim < nDim; iDim++)
     Streamwise_Coord_Vector[iDim] = config->GetPeriodic_Translation(0)[iDim];
@@ -737,6 +739,7 @@
 
   /* Value of prescribed pressure drop which results in an artificial body force vector. */
   const su2double delta_p = SPvals.Streamwise_Periodic_PressureDrop;
+  auto implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
 
   for (unsigned short iVar = 0; iVar < nVar; iVar++) residual[iVar] = 0.0;
 
@@ -749,23 +752,65 @@
   /*--- Compute the periodic temperature contribution to the energy equation, if energy equation is considered ---*/
   if (energy && streamwisePeriodic_temperature) {
 
-    scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow * DensityInc_i / (SPvals.Streamwise_Periodic_MassFlow * norm2_translation);
+    /*--- Compute scalar factors, if there is a HEAT FLUX BC ---*/
+    if (bool_heat_flux_bc) {
+      scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow * DensityInc_i / (SPvals.Streamwise_Periodic_MassFlow * norm2_translation);
 
-    /*--- Compute scalar-product dot_prod(v*t) ---*/
-    dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, &V_i[1]);
+      /*--- Compute scalar-product dot_prod(v*t) ---*/
+      dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, &V_i[1]);
+    }
+
+    /*--- Compute scalar factors, if there is an ISOTHERMAL BC ---*/
+    if (bool_isotherml_bc) {
+
+      /*--- Compute three terms associated with the source terms for iso-thermal BCs ---*/
+      // Reference Eq(20) Stalio et. al, doi:10.1115/1.2717235
+
+      // Displacement with temperature
+      dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, PrimVar_Grad_i[nDim+1]);
+
+      su2double u_i = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, &V_i[1]);
+      su2double temp_i = V_i[nDim+1];
+      su2double term1 = 2 * Thermal_Conductivity_i * dot_product;
+      su2double term2 = - SPvals.Streamwise_Periodic_LambdaL * Thermal_Conductivity_i * temp_i;
+      su2double term3 = - temp_i * u_i * DensityInc_i *  config->GetSpecific_Heat_Cp();
+      scalar_factor = SPvals.Streamwise_Periodic_LambdaL * (term1 + term2 + term3);
+      dot_product = 1.0;
+    }
 
     residual[nDim+1] = Volume * scalar_factor * dot_product;
 
-    /*--- If a RANS turbulence model is used, an additional source term, based on the eddy viscosity gradient is added. ---*/
-    if(turbulent) {
+  if (implicit) {
+
+    // /*--- Jacobian is set to zero on initialization. ---*/
 
-      /*--- Compute a scalar factor ---*/
-      scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * sqrt(norm2_translation) * Prandtl_Turb);
+    // jacobian[nDim+1][nDim+1] = Volume * scalar_factor * dot_product;
+
+  }
 
+    /*--- If a RANS turbulence model is used, an additional source term, based on the eddy viscosity gradient is added. ---*/
+    if(turbulent) {
+      if (bool_heat_flux_bc) {
+        /*--- Compute a scalar factor ---*/
+        scalar_factor = SPvals.Streamwise_Periodic_IntegratedHeatFlow / (SPvals.Streamwise_Periodic_MassFlow * sqrt(norm2_translation) * Prandtl_Turb);
+      }
+
+      if (bool_isotherml_bc) {
+        scalar_factor = (-V_i[nDim + 1] * SPvals.Streamwise_Periodic_LambdaL) * config->GetSpecific_Heat_Cp() / Prandtl_Turb;
+      }
+
       /*--- Compute scalar product between periodic translation vector and eddy viscosity gradient. ---*/
       dot_product = GeometryToolbox::DotProduct(nDim, Streamwise_Coord_Vector, AuxVar_Grad_i[0]);
 
       residual[nDim+1] -= Volume * scalar_factor * dot_product;
+
+        if (implicit) {
+
+        // /*--- Jacobian is set to zero on initialization. ---*/
+
+        // jacobian[nDim+1][nDim+1] -= Volume * scalar_factor * dot_product;
+
+        }
     } // if turbulent
   } // if energy
 

SU2_CFD/src/output/CFlowIncOutput.cpp

@@ -178,6 +178,8 @@ void CFlowIncOutput::SetHistoryOutputFields(CConfig *config){
     AddHistoryOutput("STREAMWISE_MASSFLOW", "SWMassflow", ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Massflow in streamwise periodic flow");
     AddHistoryOutput("STREAMWISE_DP",       "SWDeltaP",   ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Pressure drop in streamwise periodic flow");
     AddHistoryOutput("STREAMWISE_HEAT",     "SWHeat",     ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Integrated heat for streamwise periodic flow");
+    if (config->GetnMarker_Isothermal() > 0)
+      AddHistoryOutput("STREAMWISE_LAMBDAL", "SWLambdaL", ScreenOutputFormat::FIXED, "STREAMWISE_PERIODIC", "Exponential factor for iso-thermal temperature gradient");
   }
 
   AddAnalyzeSurfaceOutput(config);
@@ -258,10 +260,12 @@ void CFlowIncOutput::LoadHistoryData(CConfig *config, CGeometry *geometry, CSolv
 
   LoadHistoryDataScalar(config, solver);
 
-  if(streamwisePeriodic) {
+  if (streamwisePeriodic) {
     SetHistoryOutputValue("STREAMWISE_MASSFLOW", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_MassFlow);
     SetHistoryOutputValue("STREAMWISE_DP", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_PressureDrop);
     SetHistoryOutputValue("STREAMWISE_HEAT", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_IntegratedHeatFlow);
+    if (config->GetnMarker_Isothermal() > 0)
+      SetHistoryOutputValue("STREAMWISE_LAMBDAL", flow_solver->GetStreamwisePeriodicValues().Streamwise_Periodic_LambdaL);
   }
 
   /*--- Set the analyse surface history values --- */

SU2_CFD/src/output/CFlowOutput.cpp

@@ -2399,6 +2399,7 @@ void CFlowOutput::WriteMetaData(const CConfig *config){
 
     if(config->GetKind_Streamwise_Periodic() == ENUM_STREAMWISE_PERIODIC::MASSFLOW) {
       meta_file << "STREAMWISE_PERIODIC_PRESSURE_DROP=" << GetHistoryFieldValue("STREAMWISE_DP") << endl;
+      meta_file << "STREAMWISE_PERIODIC_LAMBDAL=" << GetHistoryFieldValue("STREAMWISE_LAMBDAL") << endl;
     }
   }
 

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

@@ -87,7 +87,7 @@
 
   /*--- Allocate extra solution variables, if any are in use. ---*/
 
-  const auto nVarExtra = direct_solver->RegisterSolutionExtra(true, config);
+  const auto nVarExtra = direct_solver->RegisterSolutionExtra(true, config); // TODO: Do we need this? I think no. Nitish
   nodes->AllocateAdjointSolutionExtra(nVarExtra);
 
   switch(KindDirect_Solver){
@@ -162,7 +162,7 @@
   /*--- Boolean true indicates that an input is registered ---*/
   direct_solver->GetNodes()->RegisterSolution(true);
 
-  direct_solver->RegisterSolutionExtra(true, config);
+//   direct_solver->RegisterSolutionExtra(true, config);
 
   if (time_n_needed)
     direct_solver->GetNodes()->RegisterSolution_time_n();
@@ -261,6 +261,8 @@
     config->SetIncPressureOut_BC(BPressure);
     config->SetIncTemperature_BC(Temperature);
 
+    direct_solver->RegisterSolutionExtra(reset, config);
+
   }
 
   /*--- Register incompressible radiation values as input ---*/
@@ -299,7 +301,7 @@
 
   direct_solver->GetNodes()->RegisterSolution(false);
 
-  direct_solver->RegisterSolutionExtra(false, config);
+//   direct_solver->RegisterSolutionExtra(false, config);
 }
 
 void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *config, bool CrossTerm) {