Skip to content

cleanup muscl ramp #2678

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pcarruscag merged 4 commits into from
Jan 4, 2026
+62 −71
Merged

cleanup muscl ramp #2678

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
24cc96c
cleanup
pcarruscag Jan 4, 2026
ad55dc2
fix windows build
pcarruscag Jan 4, 2026
a058ea3
fix AD issues?
pcarruscag Jan 4, 2026
8643929
no auto
pcarruscag Jan 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 14 additions & 14 deletions Common/include/option_structure.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,23 +91,23 @@ const unsigned int ZONE_0 = 0; /*!< \brief Definition of the first grid domain.
const unsigned int ZONE_1 = 1; /*!< \brief Definition of the second grid domain. */
const unsigned int INST_0 = 0; /*!< \brief Definition of the first instance per grid level. */

const su2double STANDARD_GRAVITY = 9.80665; /*!< \brief Acceleration due to gravity at surface of earth. */
const su2double UNIVERSAL_GAS_CONSTANT = 8.3144598; /*!< \brief Universal gas constant in J/(mol*K) */
const su2double BOLTZMANN_CONSTANT = 1.3806503E-23; /*!< \brief Boltzmann's constant [J K^-1] */
const su2double AVOGAD_CONSTANT = 6.0221415E26; /*!< \brief Avogadro's constant, number of particles in one kmole. */
const su2double FUND_ELEC_CHARGE_CGS = 4.8032047E-10; /*!< \brief Fundamental electric charge in CGS units, cm^(3/2) g^(1/2) s^(-1). */
constexpr passivedouble STANDARD_GRAVITY = 9.80665; /*!< \brief Acceleration due to gravity at surface of earth. */
constexpr passivedouble UNIVERSAL_GAS_CONSTANT = 8.3144598; /*!< \brief Universal gas constant in J/(mol*K) */
constexpr passivedouble BOLTZMANN_CONSTANT = 1.3806503E-23; /*!< \brief Boltzmann's constant [J K^-1] */
constexpr passivedouble AVOGAD_CONSTANT = 6.0221415E26; /*!< \brief Avogadro's constant, number of particles in one kmole. */
constexpr passivedouble FUND_ELEC_CHARGE_CGS = 4.8032047E-10; /*!< \brief Fundamental electric charge in CGS units, cm^(3/2) g^(1/2) s^(-1). */

const su2double EPS = 1.0E-16; /*!< \brief Error scale. */
const su2double TURB_EPS = 1.0E-16; /*!< \brief Turbulent Error scale. */
constexpr passivedouble EPS = 1.0E-16; /*!< \brief Error scale. */
constexpr passivedouble TURB_EPS = 1.0E-16; /*!< \brief Turbulent Error scale. */

const su2double ONE2 = 0.5; /*!< \brief One divided by two. */
const su2double ONE3 = 1.0 / 3.0; /*!< \brief One divided by three. */
const su2double TWO3 = 2.0 / 3.0; /*!< \brief Two divided by three. */
const su2double FOUR3 = 4.0 / 3.0; /*!< \brief Four divided by three. */
constexpr passivedouble ONE2 = 0.5; /*!< \brief One divided by two. */
constexpr passivedouble ONE3 = 1.0 / 3.0; /*!< \brief One divided by three. */
constexpr passivedouble TWO3 = 2.0 / 3.0; /*!< \brief Two divided by three. */
constexpr passivedouble FOUR3 = 4.0 / 3.0; /*!< \brief Four divided by three. */

const su2double PI_NUMBER = 4.0 * atan(1.0); /*!< \brief Pi number. */
constexpr passivedouble PI_NUMBER = 3.14159265358979323846; /*!< \brief Pi number (not using M_PI to avoid Windows issues). */

const su2double STEFAN_BOLTZMANN = 5.670367E-08; /*!< \brief Stefan-Boltzmann constant in W/(m^2*K^4). */
constexpr passivedouble STEFAN_BOLTZMANN = 5.670367E-08; /*!< \brief Stefan-Boltzmann constant in W/(m^2*K^4). */

const int MASTER_NODE = 0; /*!< \brief Master node for MPI parallelization. */
const int SINGLE_NODE = 1; /*!< \brief There is only a node in the MPI parallelization. */
Expand Down Expand Up @@ -195,7 +195,7 @@ const int SU2_CONN_SIZE = 10; /*!< \brief Size of the connectivity array that
that we read from a mesh file in the format [[globalID vtkType n0 n1 n2 n3 n4 n5 n6 n7 n8]. */
const int SU2_CONN_SKIP = 2; /*!< \brief Offset to skip the globalID and VTK type at the start of the element connectivity list for each CGNS element. */

const su2double COLORING_EFF_THRESH = 0.875; /*!< \brief Below this value fallback strategies are used instead. */
constexpr passivedouble COLORING_EFF_THRESH = 0.875; /*!< \brief Below this value fallback strategies are used instead. */

/*--- All temperature polynomial fits for the fluid models currently
assume a quartic form (5 coefficients). For example,
Expand Down
47 changes: 21 additions & 26 deletions SU2_CFD/include/numerics_simd/flow/convection/common.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,12 +34,12 @@

/*!
* \brief Blended difference for U-MUSCL reconstruction.
* \param[in] grad_proj - Gradient projection at point i: dot(grad_i, vector_ij).
* \param[in] gradProj - Gradient projection at point i: dot(grad_i, vector_ij).
* \param[in] delta - Centered difference: V_j - V_i.
* \param[in] kappa - Blending parameter.
* \return Blended difference for reconstruction from point i.
*/
FORCEINLINE Double umusclProjection(Double grad_proj,
FORCEINLINE Double umusclProjection(Double gradProj,
Double delta,
Double kappa) {
/*-------------------------------------------------------------------*/
Expand All @@ -51,7 +51,7 @@ FORCEINLINE Double umusclProjection(Double grad_proj,
/*--- To maintain proper scaling for edge limiters, the result of ---*/
/*--- this function is 0.5 * dV_ij^kap. ---*/
/*-------------------------------------------------------------------*/
return (1.0 - kappa) * grad_proj + kappa * delta;
return (1.0 - kappa) * gradProj + kappa * delta;
}

/*!
Expand All @@ -62,7 +62,7 @@ FORCEINLINE Double umusclProjection(Double grad_proj,
* \param[in] delta - Centered difference: V_j - V_i.
* \param[in] iVar - Variable index.
* \param[in] kappa - Blending coefficient.
* \param[in] relax - Newton-Krylov relaxation.
* \param[in] umusclRamp - MUSCL 1st-2nd order ramp times Newton-Krylov relaxation.
* \return Variable reconstructed from point i.
*/
template<class GradType, size_t nDim>
Expand All @@ -71,10 +71,9 @@ FORCEINLINE Double musclReconstruction(const GradType& grad,
const Double delta,
size_t iVar,
Double kappa,
Double relax,
Double ramp_val) {
Double umusclRamp) {
const Double proj = dot(grad[iVar], vector_ij);
return ramp_val * relax * umusclProjection(proj, delta, kappa);
return umusclRamp * umusclProjection(proj, delta, kappa);
}

/*!
Expand All @@ -87,8 +86,7 @@ FORCEINLINE void musclUnlimited(Int iPoint,
const Gradient_t& gradient,
CPair<VarType>& V,
Double kappa,
Double relax,
Double ramp_val) {
Double umusclRamp) {
constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;

auto grad_i = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
Expand All @@ -99,8 +97,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, ramp_val);
const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax, ramp_val);
const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, umusclRamp);
const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, umusclRamp);

/*--- Apply reconstruction: V_L = V_i + 0.5 * dV_ij^kap ---*/
V.i.all(iVar) += 0.5 * proj_i;
Expand All @@ -119,8 +117,7 @@ FORCEINLINE void musclPointLimited(Int iPoint,
const Gradient_t& gradient,
CPair<VarType>& V,
Double kappa,
Double relax,
Double ramp_val) {
Double umusclRamp) {
constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;

auto lim_i = gatherVariables<nVarGrad>(iPoint, limiter);
Expand All @@ -134,8 +131,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, ramp_val);
const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax, ramp_val);
const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, umusclRamp);
const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, umusclRamp);

/*--- Apply reconstruction: V_L = V_i + 0.5 * lim * dV_ij^kap ---*/
V.i.all(iVar) += 0.5 * lim_i(iVar) * proj_i;
Expand All @@ -153,8 +150,7 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
const Gradient_t& gradient,
CPair<VarType>& V,
Double kappa,
Double relax,
Double ramp_val) {
Double umusclRamp) {
constexpr auto nVarGrad = nVarGrad_ > 0 ? nVarGrad_ : VarType::nVar;

auto grad_i = gatherVariables<nVarGrad,nDim>(iPoint, gradient);
Expand All @@ -166,8 +162,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, ramp_val);
const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, relax, ramp_val);
const Double proj_i = musclReconstruction(grad_i, vector_ij, delta_ij, iVar, kappa, umusclRamp);
const Double proj_j = musclReconstruction(grad_j, vector_ij, delta_ij, iVar, kappa, umusclRamp);

/// TODO: Customize the limiter function.
const Double lim_i = (delta_ij_2 + proj_i*delta_ij) / (pow(proj_i,2) + delta_ij_2);
Expand All @@ -187,7 +183,7 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
* \param[in] gasConst - Specific gas constant.
* \param[in] muscl - If true, reconstruct, else simply fetch.
* \param[in] kappa - Blending coefficient for MUSCL reconstruction.
* \param[in] relax - Newton-Krylov relaxation.
* \param[in] umusclRamp - MUSCL 1st-2nd order ramp times Newton-Krylov relaxation.
* \param[in] limiterType - Type of flux limiter.
* \param[in] V1st - Pair of compressible flow primitives for nodes i,j.
* \param[in] vector_ij - Distance vector from i to j.
Expand All @@ -200,12 +196,11 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
const su2double& gasConst,
bool muscl,
const su2double& kappa,
const su2double& relax,
const su2double& umusclRamp,
LIMITER limiterType,
const CPair<PrimVarType>& V1st,
const VectorDbl<nDim>& vector_ij,
const VariableType& solution,
const su2double& ramp_val) {
const VariableType& solution) {
static_assert(ReconVarType::nVar <= PrimVarType::nVar);

const auto& gradients = solution.GetGradient_Reconstruction();
Expand All @@ -223,13 +218,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, ramp_val);
musclUnlimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, kappa, umusclRamp);
break;
case LIMITER::VAN_ALBADA_EDGE:
musclEdgeLimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, kappa, relax, ramp_val);
musclEdgeLimited<nVarGrad>(iPoint, jPoint, vector_ij, gradients, V, kappa, umusclRamp);
break;
default:
musclPointLimited<nVarGrad>(iPoint, jPoint, vector_ij, limiters, gradients, V, kappa, relax, ramp_val);
musclPointLimited<nVarGrad>(iPoint, jPoint, vector_ij, limiters, gradients, V, kappa, umusclRamp);
break;
}
/*--- Recompute density using the reconstructed pressure and temperature. ---*/
Expand Down
4 changes: 2 additions & 2 deletions SU2_CFD/include/numerics_simd/flow/convection/roe.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ class CRoeBase : public Base {
AD::StartPreacc();

const bool implicit = (config.GetKind_TimeIntScheme() == EULER_IMPLICIT);
const auto nkRelax = config.GetNewtonKrylovRelaxation();
const su2double umusclRamp = config.GetMUSCLRampValue() * config.GetNewtonKrylovRelaxation();
const auto& solution = static_cast<const CEulerVariable&>(solution_);

const auto iPoint = geometry.edges->GetNode(iEdge,0);
Expand All @@ -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, config.GetMUSCLRampValue());
iEdge, iPoint, jPoint, gamma, gasConst, muscl, umusclKappa, umusclRamp, typeLimiter, V1st, vector_ij, solution);

/*--- Compute conservative variables. ---*/

Expand Down
9 changes: 5 additions & 4 deletions SU2_CFD/include/solvers/CScalarSolver.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,7 @@ void CScalarSolver<VariableType>::Upwind_Residual(CGeometry* geometry, CSolver**
/*--- U-MUSCL reconstruction ---*/
const su2double kappa = config->GetMUSCL_Kappa();
const su2double kappaFlow = config->GetMUSCL_Kappa_Flow();
const su2double musclRamp = config->GetMUSCLRampValue();

auto* flowNodes = su2staticcast_p<CFlowVariable*>(solver_container[FLOW_SOL]->GetNodes());
const auto& edgeMassFluxes = *(solver_container[FLOW_SOL]->GetEdgeMassFluxes());
Expand Down Expand Up @@ -222,8 +223,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, config->GetMUSCLRampValue());
su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappaFlow, config->GetMUSCLRampValue());
su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappaFlow, musclRamp);
su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappaFlow, musclRamp);

if (limiterFlow) {
Project_Grad_i *= Limiter_i[iVar];
Expand Down Expand Up @@ -251,8 +252,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, config->GetMUSCLRampValue());
su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, U_ij, kappa, config->GetMUSCLRampValue());
su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, U_ij, kappa, musclRamp);
su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, U_ij, kappa, musclRamp);

if (limiter) {
Project_Grad_i *= Limiter_i[iVar];
Expand Down
5 changes: 3 additions & 2 deletions SU2_CFD/include/solvers/CSolver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -587,10 +587,11 @@ 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
* \param[in] ramp_val - Value of the 1st-2nd order MUSCL ramp.
* \return - Projected variable.
*/
inline su2double MUSCL_Reconstruction(const su2double* grad, const su2double* vector_ij, su2double delta_ij, su2double kappa, su2double ramp_val) {
FORCEINLINE su2double MUSCL_Reconstruction(const su2double* grad, const su2double* vector_ij, su2double delta_ij,
su2double kappa, su2double ramp_val) const {
su2double project_grad = GeometryToolbox::DotProduct(nDim, grad, vector_ij);
return ramp_val * LimiterHelpers<>::umusclProjection(project_grad, delta_ij, kappa);
}
Expand Down
12 changes: 4 additions & 8 deletions SU2_CFD/src/integration/CNewtonIntegration.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,15 +184,11 @@ void CNewtonIntegration::MultiGrid_Iteration(CGeometry ****geometry_, CSolver **

if (!setup) { Setup(); setup = true; }

// Ramp from 1st to 2nd order during the startup.
su2double baseNkRelaxation = 1;
if (startupPeriod && startupIters > 0 && !config->GetRestart()) {
baseNkRelaxation = su2double(startupIters - iter) / startupIters;
}
config->SetNewtonKrylovRelaxation(baseNkRelaxation);
/*--- Remove NK relaxation to compute the current residual. ---*/
config->SetNewtonKrylovRelaxation(1.0);

// When using NK relaxation (not fully 2nd order Jacobian products) we need an additional
// residual evaluation that is used as the reference for finite differences.
/*--- When using NK relaxation (not fully 2nd order Jacobian products) we need an additional
* residual evaluation that is used as the reference for finite differences. ---*/
LinSysRes0 = (!startupPeriod && nkRelaxation < 1) ? &LinSysResRelax : &LinSysRes;

SU2_OMP_PARALLEL_(if(solvers[FLOW_SOL]->GetHasHybridParallel())) {
Expand Down
2 changes: 1 addition & 1 deletion SU2_CFD/src/iteration/CFluidIteration.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -348,7 +348,7 @@ void CFluidIteration::UpdateRamp(CGeometry**** geometry_container, CConfig** con
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));
config->SetMUSCLRampValue(std::pow((0.5 * (1 - cos(PI_NUMBER * std::min(1.0, iterFrac)))), RampMUSCLParam.RampMUSCLPower));
break;
default:
break;
Expand Down
6 changes: 3 additions & 3 deletions SU2_CFD/src/solvers/CEulerSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1799,7 +1799,6 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
}

const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
const su2double nkRelax = config->GetNewtonKrylovRelaxation();

const bool roe_turkel = (config->GetKind_Upwind_Flow() == UPWIND::TURKEL);
const auto kind_dissipation = config->GetKind_RoeLowDiss();
Expand All @@ -1809,6 +1808,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE);

const su2double kappa = config->GetMUSCL_Kappa_Flow();
const su2double musclRamp = config->GetMUSCLRampValue() * config->GetNewtonKrylovRelaxation();

/*--- Non-physical counter. ---*/
unsigned long counter_local = 0;
Expand Down Expand Up @@ -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, config->GetMUSCLRampValue());
const su2double Project_Grad_j = nkRelax * MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappa, config->GetMUSCLRampValue());
const su2double Project_Grad_i = MUSCL_Reconstruction(Gradient_i[iVar], Vector_ij, V_ij, kappa, musclRamp);
const su2double Project_Grad_j = MUSCL_Reconstruction(Gradient_j[iVar], Vector_ij, V_ij, kappa, musclRamp);

su2double lim_i = 1.0;
su2double lim_j = 1.0;
Expand Down
Loading
Loading