Merge pull request #2639 from joshkellyjak/feature_MUSCL_ramp

MUSCL Ramp options

Author: pcarruscag

Common/include/CConfig.hpp

@@ -1106,6 +1106,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 */
@@ -1122,6 +1123,13 @@ 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. */
 
+  struct CMUSCLRampParam {
+    su2double RampMUSCLPower; /*!< \brief Exponent by which to raise the MUSCL ramp to the power of */
+    MUSCL_RAMP_TYPE Kind_MUSCLRamp; /*!< \brief The kind of MUSCL ramp */
+    unsigned long rampMUSCLCoeff[3];     /*!< \brief ramp MUSCL value coefficients for the COption class. */
+  } RampMUSCLParam;
+  su2double rampMUSCLValue; /*!< \brief Current value of the MUSCL ramp */
+
   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. */
@@ -1352,6 +1360,8 @@ class CConfig {
 
   void addUShortArrayOption(const string& name, int size, bool allow_fewer, unsigned short* option_field);
 
+  void addULongArrayOption(const string& name, int size, bool allow_fewer, unsigned long* option_field);
+
   void addDoubleListOption(const string& name, unsigned short & size, su2double * & option_field);
 
   void addShortListOption(const string& name, unsigned short & size, short * & option_field);
@@ -5208,6 +5218,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
@@ -5218,6 +5234,23 @@ class CConfig {
     else return 0;
   };
 
+  /*!
+   * \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 paramaters.
+   * \return Ramp MUSCL kind
+  */
+  const CMUSCLRampParam& GetMUSCLRampParam(void) const { return RampMUSCLParam; }
+
   /*!
    * \brief Generic interface for setting monitor outlet values for the ramp.
    */

Common/include/option_structure.hpp

@@ -1959,7 +1959,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

@@ -409,6 +409,13 @@ void CConfig::addULongListOption(const string& name, unsigned short & size, unsi
   option_map.insert(pair<string, COptionBase *>(name, val));
 }
 
+void CConfig::addULongArrayOption(const string& name, const int size, const bool allow_fewer, unsigned long* option_field) {
+  assert(option_map.find(name) == option_map.end());
+  all_options.insert(pair<string, bool>(name, true));
+  COptionBase* val = new COptionArray<unsigned long>(name, size, allow_fewer, option_field);
+  option_map.insert(pair<string, COptionBase *>(name, val));
+}
+
 void CConfig::addStringListOption(const string& name, unsigned short & num_marker, string* & option_field) {
   assert(option_map.find(name) == option_map.end());
   all_options.insert(pair<string, bool>(name, true));
@@ -2009,6 +2016,16 @@ 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 */
+  RampMUSCLParam.rampMUSCLCoeff[0] = 0.0; RampMUSCLParam.rampMUSCLCoeff[1] = 1.0; RampMUSCLParam.rampMUSCLCoeff[2] = 500.0;
+  addULongArrayOption("RAMP_MUSCL_COEFF", 3, false, RampMUSCLParam.rampMUSCLCoeff);
+  /*!\brief RAMP_MUSCL_POWER \n DESRCIPTION: Exponent of the MUSCL ramp formulation */
+  addDoubleOption("RAMP_MUSCL_POWER", RampMUSCLParam.RampMUSCLPower, 1.0);
+  /*!\brief KIND_MUSCL_RAMP \n DESCRIPTION: The kind of MUSCL Ramp to be applied */
+  addEnumOption("KIND_MUSCL_RAMP", RampMUSCLParam.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);
@@ -4502,6 +4519,13 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     }
   }
 
+  if(RampMUSCL && !DiscreteAdjoint){
+    if (RampMUSCLParam.RampMUSCLPower <= 0.0) SU2_MPI::Error("RAMP_MUSCL_POWER cannot be less than or equal to zero!", CURRENT_FUNCTION);
+    rampMUSCLValue = 0.0;
+  } else {
+    rampMUSCLValue = 1.0;
+  }
+
   /*--- Check on extra Relaxation factor for Giles---*/
   if(extrarelfac[1] > 0.5){
     extrarelfac[1] = 0.5;
@@ -6993,6 +7017,21 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
 
     auto PrintLimiterInfo = [&](const LIMITER kind_limiter, const su2double kappa) {
       cout << "Second order integration in space, with slope limiter.\n";
+      if (RampMUSCL) {
+        cout << "Ramping MUSCL sheme from first to second order starting at iter " << RampMUSCLParam.rampMUSCLCoeff[RAMP_COEFF::INITIAL_VALUE]
+          << ", ending at iter " << RampMUSCLParam.rampMUSCLCoeff[RAMP_COEFF::FINAL_ITER] + RampMUSCLParam.rampMUSCLCoeff[RAMP_COEFF::INITIAL_VALUE]
+          << ", updating every " << RampMUSCLParam.rampMUSCLCoeff[RAMP_COEFF::UPDATE_FREQ] << " iterations." << endl;
+        string MUSCLRampType;
+        switch (RampMUSCLParam.Kind_MUSCLRamp) {
+          case MUSCL_RAMP_TYPE::ITERATION:
+            MUSCLRampType = "linear";
+            break;
+          case MUSCL_RAMP_TYPE::SMOOTH_FUNCTION:
+            MUSCLRampType = "cosine";
+            break;
+        }
+        cout << "Ramp applied according to a " << MUSCLRampType << " function, raised to the power " << RampMUSCLParam.RampMUSCLPower << "." << endl;
+      }
       if (kappa != 0.0) cout << "U-MUSCL reconstruction, with coefficient: " << kappa << ".\n";
       switch (kind_limiter) {
         case LIMITER::NONE:

SU2_CFD/include/numerics_simd/flow/convection/common.hpp

@@ -71,9 +71,10 @@ FORCEINLINE Double musclReconstruction(const GradType& grad,
                                        const Double delta,
                                        size_t iVar,
                                        Double kappa,
-                                       Double relax) {
+                                       Double relax,
+                                       Double ramp_val) {
   const Double proj = dot(grad[iVar], vector_ij);
-  return relax * umusclProjection(proj, delta, kappa);
+  return ramp_val * relax * umusclProjection(proj, delta, kappa);
 }
 
 /*!
@@ -86,7 +87,8 @@ FORCEINLINE void musclUnlimited(Int iPoint,
                                 const Gradient_t& gradient,
                                 CPair<VarType>& V,
                                 Double kappa,
-                                Double relax) {
+                                Double relax,
+                                Double ramp_val) {
   constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
 
   auto grad_i = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
@@ -97,8 +99,8 @@ FORCEINLINE void musclUnlimited(Int iPoint,
     const Double delta_ij = V.j.all(iVar) - V.i.all(iVar);
 
     /*--- U-MUSCL reconstructed variables ---*/
-    const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, relax);
-    const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax);
+    const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, relax, ramp_val);
+    const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax, ramp_val);
 
     /*--- Apply reconstruction: V_L = V_i + 0.5 * dV_ij^kap ---*/
     V.i.all(iVar) += 0.5 * proj_i;
@@ -117,7 +119,8 @@ FORCEINLINE void musclPointLimited(Int iPoint,
                                    const Gradient_t& gradient,
                                    CPair<VarType>& V,
                                    Double kappa,
-                                   Double relax) {
+                                   Double relax,
+                                   Double ramp_val) {
   constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
 
   auto lim_i = gatherVariables<nVarGrad>(iPoint, limiter);
@@ -131,8 +134,8 @@ FORCEINLINE void musclPointLimited(Int iPoint,
     const Double delta_ij = V.j.all(iVar) - V.i.all(iVar);
 
     /*--- U-MUSCL reconstructed variables ---*/
-    const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, relax);
-    const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax);
+    const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, relax, ramp_val);
+    const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax, ramp_val);
 
     /*--- Apply reconstruction: V_L = V_i + 0.5 * lim * dV_ij^kap ---*/
     V.i.all(iVar) += 0.5 * lim_i(iVar) * proj_i;
@@ -150,7 +153,8 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
                                   const Gradient_t& gradient,
                                   CPair<VarType>& V,
                                   Double kappa,
-                                  Double relax) {
+                                  Double relax,
+                                  Double ramp_val) {
   constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;
 
   auto grad_i = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
@@ -162,8 +166,8 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
     const Double delta_ij_2 = pow(delta_ij, 2) + 1e-6;
 
     /*--- U-MUSCL reconstructed variables ---*/
-    const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, relax);
-    const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax);
+    const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, relax, ramp_val);
+    const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax, ramp_val);
 
     /// TODO: Customize the limiter function.
     const Double lim_i = (delta_ij_2 + proj_i*delta_ij) / (pow(proj_i,2) + delta_ij_2);
@@ -200,7 +204,8 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
                                                       LIMITER limiterType,
                                                       const CPair<PrimVarType>& V1st,
                                                       const VectorDbl<nDim>& vector_ij,
-                                                      const VariableType& solution) {
+                                                      const VariableType& solution,
+                                                      const su2double& ramp_val) {
   static_assert(ReconVarType::nVar <= PrimVarType::nVar);
 
   const auto& gradients = solution.GetGradient_Reconstruction();
@@ -218,13 +223,13 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
     constexpr auto nVarGrad = ReconVarType::nVar - 2;
     switch (limiterType) {
     case LIMITER::NONE:
-      musclUnlimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, kappa, relax);
+      musclUnlimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, kappa, relax, ramp_val);
       break;
     case LIMITER::VAN_ALBADA_EDGE:
-      musclEdgeLimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, kappa, relax);
+      musclEdgeLimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, kappa, relax, ramp_val);
       break;
     default:
-      musclPointLimited<nVarGrad>(iPoint, jPoint, vector_ij, limiters, gradients, V, kappa, relax);
+      musclPointLimited<nVarGrad>(iPoint, jPoint, vector_ij, limiters, gradients, V, kappa, relax, ramp_val);
       break;
     }
     /*--- Recompute density using the reconstructed pressure and temperature. ---*/

SU2_CFD/include/numerics_simd/flow/convection/roe.hpp

@@ -123,7 +123,7 @@ class CRoeBase : public Base {
 
     /*--- Recompute density and enthalpy instead of reconstructing. ---*/
     auto V = reconstructPrimitives<CCompressiblePrimitives<nDim,nPrimVarGrad> >(
-        iEdge, iPoint, jPoint, gamma, gasConst, muscl, umusclKappa, nkRelax, typeLimiter, V1st, vector_ij, solution);
+        iEdge, iPoint, jPoint, gamma, gasConst, muscl, umusclKappa, nkRelax, typeLimiter, V1st, vector_ij, solution, config.GetMUSCLRampValue());
 
     /*--- Compute conservative variables. ---*/
 

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 = MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappaFlow, config->GetMUSCLRampValue());
+            su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappaFlow, config->GetMUSCLRampValue());
 
             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 = MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, U_ij, kappa, config->GetMUSCLRampValue());
+            su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, U_ij, kappa, config->GetMUSCLRampValue());
 
             if (limiter) {
               Project_Grad_i *= Limiter_i[iVar];

SU2_CFD/include/solvers/CSolver.hpp

@@ -587,11 +587,12 @@ 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) {
+  inline su2double MUSCL_Reconstruction(const su2double* grad, const su2double* vector_ij, su2double delta_ij, su2double kappa, su2double ramp_val) {
     su2double project_grad = GeometryToolbox::DotProduct(nDim, grad, vector_ij);
-    return LimiterHelpers<>::umusclProjection(project_grad, delta_ij, kappa);
+    return ramp_val * LimiterHelpers<>::umusclProjection(project_grad, delta_ij, 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 auto RampMUSCLParam = config->GetMUSCLRampParam();
+    const long unsigned startIter = RampMUSCLParam.rampMUSCLCoeff[RAMP_COEFF::INITIAL_VALUE];
+    const long unsigned updateFreq = RampMUSCLParam.rampMUSCLCoeff[RAMP_COEFF::UPDATE_FREQ];
+    const long unsigned rampLength = RampMUSCLParam.rampMUSCLCoeff[RAMP_COEFF::FINAL_ITER];
+    auto iterFrac = (static_cast<double>(iter - startIter)/static_cast<double>((rampLength + startIter) - startIter));
+    if (iter < startIter) return;
+    if ((iter == startIter) && (rank == MASTER_NODE)) cout << "Beginning to ramp MUSCL scheme..." << endl;
+    if ((iter % updateFreq == 0 && iter < (rampLength + startIter)) || (iter == (rampLength + startIter))) {
+      switch (RampMUSCLParam.Kind_MUSCLRamp) {
+        case MUSCL_RAMP_TYPE::ITERATION:
+          config->SetMUSCLRampValue(std::pow(std::min<double>(1.0, iterFrac), RampMUSCLParam.RampMUSCLPower));
+          break;
+        case MUSCL_RAMP_TYPE::SMOOTH_FUNCTION:
+          config->SetMUSCLRampValue(std::pow((0.5 * (1 -  cos(M_PI * std::min(1.0, iterFrac)))), RampMUSCLParam.RampMUSCLPower));
+          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

@@ -1886,8 +1886,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 = nkRelax * MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappa, config->GetMUSCLRampValue());
+        const su2double Project_Grad_j = nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappa, config->GetMUSCLRampValue());
 
         su2double lim_i = 1.0;
         su2double lim_j = 1.0;