Added MUSCL Ramps options

Author: joshkellyjak

Common/include/CConfig.hpp

@@ -1097,6 +1097,7 @@ class CConfig {
   bool SpatialFourier;              /*!< \brief option for computing the fourier transforms for subsonic non-reflecting BC. */
   bool RampMotionFrame;             /*!< \brief option for ramping up or down the motion Frame values */
   bool RampOutlet;                  /*!< \brief option for ramping up or down the outlet values */
+  bool RampMUSCL;
   bool RampRotatingFrame;           /*!< \brief option for ramping up or down the motion Frame values */
   bool RampTranslationFrame;        /*!< \brief option for ramping up or down the outlet values */
   bool RampOutletMassFlow;          /*!< \brief option for ramping up or down the motion Frame values */
@@ -1113,6 +1114,11 @@ class CConfig {
   array<su2double, N_POLY_COEFFS> kt_polycoeffs{{0.0}};  /*!< \brief Array for thermal conductivity polynomial coefficients. */
   bool Body_Force;                      /*!< \brief Flag to know if a body force is included in the formulation. */
 
+  su2double rampMUSCLValue; /*!< \brief Current value of the MUSCL ramp */
+  su2double RampMUSCLPower; /*!< \brief Exponent by which to raise the MUSCL ramp to the power of */
+  MUSCL_RAMP_TYPE Kind_MUSCLRamp;
+  unsigned long *rampMUSCLCoeff;     /*!< \brief ramp MUSCL value coefficients for the COption class. */
+
   ENUM_STREAMWISE_PERIODIC Kind_Streamwise_Periodic; /*!< \brief Kind of Streamwise periodic flow (pressure drop or massflow) */
   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. */
@@ -5179,6 +5185,12 @@ class CConfig {
    */
   bool GetRampOutflow(void) const { return RampOutlet; }
 
+  /*!
+   * \brief Get MUSCL ramp option.
+   * \return Ramp MUSCL option
+  */
+  bool GetMUSCLRamp(void) const { return RampMUSCL; }
+
   /*!
    * \brief General interface for accessing ramp coefficient information
    * \return coeff for ramps
@@ -5189,6 +5201,37 @@ class CConfig {
     else return 0;
   };
 
+  /*!
+   * \brief Interface for accessing MUSCL ramp coefficient information
+   * \return coeff for ramps
+  */
+  unsigned long GetMUSCLRampCoeff(RAMP_COEFF val_coeff) const {
+    return rampMUSCLCoeff[val_coeff];
+  }
+
+  /*!
+   * \brief Set MUSCL ramp value.
+  */
+  void SetMUSCLRampValue(su2double ramp_value) { rampMUSCLValue = ramp_value; }
+
+  /*!
+   * \brief Get MUSCL ramp value.
+   * \return Ramp MUSCL value
+  */
+  su2double GetMUSCLRampValue(void) const { return rampMUSCLValue; }
+
+  /*!
+   * \brief Get MUSCL ramp power.
+   * \return Ramp MUSCL power
+  */
+  su2double GetMUSCLRampPower(void) const { return RampMUSCLPower; }
+
+  /*!
+   * \brief Get MUSCL ramp kind.
+   * \return Ramp MUSCL kind
+  */
+  MUSCL_RAMP_TYPE GetKind_MUSCLRamp(void) const { return Kind_MUSCLRamp; }
+
   /*!
    * \brief Generic interface for setting monitor outlet values for the ramp.
    */

Common/include/option_structure.hpp

@@ -1931,7 +1931,18 @@ enum TURBO_MARKER_TYPE{
 
 enum class RAMP_TYPE{
   GRID,       /*!< \brief flag for rotational/translational ramps */
-  BOUNDARY    /*!< \brief flag for pressure/mass flow ramps*/
+  BOUNDARY,   /*!< \brief flag for pressure/mass flow ramps*/
+  MUSCL       /*!< \brief flag for MUSCL ramps */
+};
+
+enum class MUSCL_RAMP_TYPE{
+  ITERATION,  /*!< \brief flag for linear iteration-based ramp */
+  SMOOTH_FUNCTION      /*!< \brief flag for smooth cosine ramp */
+};
+
+static const MapType<std::string, MUSCL_RAMP_TYPE> MUSCLRamp_Map = {
+  MakePair("ITERATION", MUSCL_RAMP_TYPE::ITERATION)
+  MakePair("SMOOTH_FUNCTION", MUSCL_RAMP_TYPE::SMOOTH_FUNCTION)
 };
 
 /*!

Common/src/CConfig.cpp

@@ -1986,6 +1986,18 @@ void CConfig::SetConfig_Options() {
   addBoolOption("MUSCL_FLOW", MUSCL_Flow, true);
   /*!\brief MUSCL_KAPPA_FLOW \n DESCRIPTION: Blending coefficient for the U-MUSCL scheme \ingroup Config*/
   addDoubleOption("MUSCL_KAPPA_FLOW", MUSCL_Kappa_Flow, 0.0);
+  /*!\brief RAMP_MUSCL \n DESCRIPTION: Enable ramping of the MUSCL scheme from 1st to 2nd order using specified method*/
+  addBoolOption("RAMP_MUSCL", RampMUSCL, false);
+  /*! brief RAMP_OUTLET_COEFF \n DESCRIPTION: the 1st coeff is the ramp start iteration,
+   * the 2nd coeff is the iteration update frequenct, 3rd coeff is the total number of iterations */
+  unsigned short nMUSCLCoeffs = 3;
+  rampMUSCLCoeff = new unsigned long [nMUSCLCoeffs];
+  rampMUSCLCoeff[0] = 0.0; rampMUSCLCoeff[1] = 1.0; rampMUSCLCoeff[2] = 500.0;
+  addULongListOption("RAMP_MUSCL_COEFF", nMUSCLCoeffs, rampMUSCLCoeff);
+  /*!\brief RAMP_MUSCL_POWER \n DESRCIPTION: Exponent of the MUSCL ramp formulation */
+  addDoubleOption("RAMP_MUSCL_POWER", RampMUSCLPower, 1.0);
+  /*!\brief KIND_MUSCL_RAMP \n DESCRIPTION: The kind of MUSCL Ramp to be applied */
+  addEnumOption("KIND_MUSCL_RAMP", Kind_MUSCLRamp, MUSCLRamp_Map, MUSCL_RAMP_TYPE::ITERATION);
   /*!\brief SLOPE_LIMITER_FLOW
    * DESCRIPTION: Slope limiter for the direct solution. \n OPTIONS: See \link Limiter_Map \endlink \n DEFAULT VENKATAKRISHNAN \ingroup Config*/
   addEnumOption("SLOPE_LIMITER_FLOW", Kind_SlopeLimit_Flow, Limiter_Map, LIMITER::VENKATAKRISHNAN);
@@ -4479,6 +4491,12 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     }
   }
 
+  if(RampMUSCL && !DiscreteAdjoint){
+    rampMUSCLValue = 0.0;
+  } else {
+    rampMUSCLValue = 1.0;
+  }
+
   /*--- Check on extra Relaxation factor for Giles---*/
   if(extrarelfac[1] > 0.5){
     extrarelfac[1] = 0.5;

SU2_CFD/include/solvers/CScalarSolver.inl

@@ -222,8 +222,8 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
           for (auto iVar = 0u; iVar < solver_container[FLOW_SOL]->GetnPrimVarGrad(); iVar++) {
             const su2double V_ij = V_j[iVar] - V_i[iVar];
 
-            su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappaFlow);
-            su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappaFlow);
+            su2double Project_Grad_i = config->GetMUSCLRampValue() * MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappaFlow);
+            su2double Project_Grad_j = config->GetMUSCLRampValue() * MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappaFlow);
 
             if (limiterFlow) {
               Project_Grad_i *= Limiter_i[iVar];
@@ -251,8 +251,8 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
           for (auto iVar = 0u; iVar < nVar; iVar++) {
             const su2double U_ij = Scalar_j[iVar] - Scalar_i[iVar];
 
-            su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, U_ij, kappa);
-            su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, U_ij, kappa);
+            su2double Project_Grad_i = config->GetMUSCLRampValue() * MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, U_ij, kappa);
+            su2double Project_Grad_j = config->GetMUSCLRampValue() * MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, U_ij, kappa);
 
             if (limiter) {
               Project_Grad_i *= Limiter_i[iVar];

SU2_CFD/include/solvers/CSolver.hpp

@@ -587,6 +587,7 @@ class CSolver {
    * \param[in] vector_ij - Distance vector.
    * \param[in] delta_ij - Centered difference.
    * \param[in] kappa - Blending coefficient for U-MUSCL reconstruction.
+   * \param[in] ramp_val - Value of the ramp
    * \return - Projected variable.
    */
   inline su2double MUSCL_Reconstruction(const su2double* grad, const su2double* vector_ij, su2double delta_ij, su2double kappa) {

SU2_CFD/src/iteration/CFluidIteration.cpp

@@ -241,6 +241,9 @@ bool CFluidIteration::Monitor(COutput* output, CIntegration**** integration, CGe
   if (config[val_iZone]->GetRampOutflow())
       UpdateRamp(geometry, config, config[val_iZone]->GetInnerIter(), val_iZone, RAMP_TYPE::BOUNDARY);
 
+  if (config[val_iZone]->GetMUSCLRamp())
+    UpdateRamp(geometry, config, config[val_iZone]->GetInnerIter(), val_iZone, RAMP_TYPE::MUSCL);
+
   output->SetHistoryOutput(geometry[val_iZone][val_iInst][MESH_0], solver[val_iZone][val_iInst][MESH_0],
                            config[val_iZone], config[val_iZone]->GetTimeIter(), config[val_iZone]->GetOuterIter(),
                            config[val_iZone]->GetInnerIter());
@@ -330,6 +333,29 @@ void CFluidIteration::UpdateRamp(CGeometry**** geometry_container, CConfig** con
       }
     }
   }
+
+  if (ramp_flag == RAMP_TYPE::MUSCL) {
+    const long unsigned startIter = config->GetMUSCLRampCoeff(RAMP_COEFF::INITIAL_VALUE);
+    const long unsigned updateFreq = config->GetMUSCLRampCoeff(RAMP_COEFF::UPDATE_FREQ);
+    const long unsigned finalIter = config->GetMUSCLRampCoeff(RAMP_COEFF::FINAL_ITER);
+    const auto power = config->GetMUSCLRampPower();
+    auto iterFrac = (static_cast<double>(iter - startIter)/static_cast<double>(finalIter - startIter));
+    if (iter < startIter) return;
+    if ((iter == startIter) && (rank == MASTER_NODE)) cout << "Beginning to ramp MUSCL scheme..." << endl;
+    if ((iter % updateFreq == 0 && iter < finalIter) || (iter == finalIter)) {
+      switch (config->GetKind_MUSCLRamp()) {
+        case MUSCL_RAMP_TYPE::ITERATION:
+          config->SetMUSCLRampValue(std::pow(std::min<double>(1.0, iterFrac), power));
+          break;
+        case MUSCL_RAMP_TYPE::SMOOTH_FUNCTION:
+          config->SetMUSCLRampValue(std::pow((0.5 * (1 -  cos(M_PI * std::min(1.0, iterFrac)))), power));
+          break;
+        default:
+          break;
+      }
+      if (rank == MASTER_NODE) cout << "MUSCL Ramp value updated. New Value: " << config->GetMUSCLRampValue() << endl;
+    }
+  }
 }
 
 void CFluidIteration::ComputeTurboPerformance(CSolver***** solver, CGeometry**** geometry_container, CConfig** config_container, unsigned long ExtIter) {

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -1881,8 +1881,8 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
       for (auto iVar = 0u; iVar < nPrimVarGrad; iVar++) {
         const su2double V_ij = V_j[iVar] - V_i[iVar];
 
-        const su2double Project_Grad_i = nkRelax * MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappa);
-        const su2double Project_Grad_j = nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappa);
+        const su2double Project_Grad_i = config->GetMUSCLRampValue() * nkRelax * MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappa);
+        const su2double Project_Grad_j = config->GetMUSCLRampValue() * nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappa);
 
         su2double lim_i = 1.0;
         su2double lim_j = 1.0;

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -1279,8 +1279,8 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
       for (auto iVar = 0u; iVar < nPrimVarGrad; iVar++) {
         const su2double V_ij = V_j[iVar] - V_i[iVar];
 
-        const su2double Project_Grad_i = nkRelax * MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappa);
-        const su2double Project_Grad_j = nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappa);
+        const su2double Project_Grad_i = config->GetMUSCLRampValue() * nkRelax * MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappa);
+        const su2double Project_Grad_j = config->GetMUSCLRampValue() * nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappa);
 
         su2double lim_i = 1.0;
         su2double lim_j = 1.0;

SU2_CFD/src/solvers/CNEMOEulerSolver.cpp

@@ -542,8 +542,8 @@ void CNEMOEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_con
       for (auto iVar = 0ul; iVar < nPrimVarGrad; iVar++) {
         const su2double V_ij = V_j[iVar] - V_i[iVar];
 
-        Project_Grad_i[iVar] = nkRelax * MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappa);
-        Project_Grad_j[iVar] = nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappa);
+        Project_Grad_i[iVar] = config->GetMUSCLRampValue() * nkRelax * MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappa);
+        Project_Grad_j[iVar] = config->GetMUSCLRampValue() * nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappa);
 
         if (limiter) {
           if (van_albada) {