Skip to content

Commit 4cd4784

Browse files
committed
introduce LSQ slope computation
1 parent 6cfdbf0 commit 4cd4784

File tree

1 file changed

+19
-17
lines changed

1 file changed

+19
-17
lines changed

SU2_CFD/src/solvers/CSolver.cpp

Lines changed: 19 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1908,19 +1908,26 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
19081908

19091909
auto slopeLSQ = [](const deque<pair<unsigned long, su2double>>& pts) {
19101910
if (pts.size() < 2) return 0.0;
1911-
if (pts.size() == 2) {
1912-
/* Two points: simple slope. */
1913-
const su2double dt = su2double(pts[1].first - pts[0].first);
1914-
if (dt < 1e-12) return 0.0; /* Avoid division by zero */
1915-
return (pts[1].second - pts[0].second) / dt;
1911+
1912+
/* Compute least-squares linear regression slope:
1913+
slope = [n*Σ(xy) - Σ(x)*Σ(y)] / [n*Σ(x²) - (Σ(x))²] */
1914+
su2double n = su2double(pts.size());
1915+
su2double sumX = 0.0, sumY = 0.0, sumXY = 0.0, sumX2 = 0.0;
1916+
1917+
for (const auto& pt : pts) {
1918+
su2double x = su2double(pt.first); // iteration number
1919+
su2double y = pt.second; // log10 residual value
1920+
sumX += x;
1921+
sumY += y;
1922+
sumXY += x * y;
1923+
sumX2 += x * x;
19161924
}
1917-
/* Three points: average of two consecutive slopes. */
1918-
const su2double dt1 = su2double(pts[1].first - pts[0].first);
1919-
const su2double dt2 = su2double(pts[2].first - pts[1].first);
1920-
if (dt1 < 1e-12 || dt2 < 1e-12) return 0.0; /* Avoid division by zero */
1921-
const su2double slope1 = (pts[1].second - pts[0].second) / dt1;
1922-
const su2double slope2 = (pts[2].second - pts[1].second) / dt2;
1923-
return 0.5 * (slope1 + slope2);
1925+
1926+
su2double denominator = n * sumX2 - sumX * sumX;
1927+
if (fabs(denominator) < 1e-12) return 0.0; // Avoid division by zero
1928+
1929+
su2double slope = (n * sumXY - sumX * sumY) / denominator;
1930+
return slope;
19241931
};
19251932

19261933
bool peaksStagnant = false, valleysStagnant = false;
@@ -2199,11 +2206,6 @@ void CSolver::AdaptCFLNumber(CGeometry **geometry,
21992206

22002207
}
22012208

2202-
//for (unsigned short iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) {
2203-
// const su2double CFLRatio = config->GetCFL(iMesh)/config->GetCFL(iMesh-1);
2204-
// MGFactor[iMesh] = MGFactor[iMesh-1]*CFLRatio;
2205-
// cout << " CFL Level " << iMesh << " initial: " << config->GetCFL(iMesh)<< endl;
2206-
// cout << " CFL Level " << iMesh-1 << " initial: " << config->GetCFL(iMesh-1) << endl;
22072209

22082210
}
22092211

0 commit comments

Comments
 (0)