From 92033815bad47f21fa3ae47cbc0a865dd852ac82 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Corentin=20Thom=C3=A9e?= <corentin.thomee@student.uliege.be>
Date: Wed, 13 Mar 2024 14:11:44 +0100
Subject: [PATCH] Handled degenrate case in Kutta

---
 hspm/src/geometry.cpp | 2 +-
 hspm/src/solver.cpp   | 5 ++++-
 python/main.py        | 9 +--------
 3 files changed, 6 insertions(+), 10 deletions(-)

diff --git a/hspm/src/geometry.cpp b/hspm/src/geometry.cpp
index e8ec926..efd5b93 100644
--- a/hspm/src/geometry.cpp
+++ b/hspm/src/geometry.cpp
@@ -30,7 +30,7 @@ void HSPM::generateNaca4DigitCoordinates(double camber, double camberPos, double
 		} else {
 			x_distr(i) = chord/2 * (cos(angle)+1);
 		}*/
-		x_distr(i) = chord * pow((1 + cos((angle + M_PI) / 2)),2);
+		x_distr(i) = chord * pow((1 + cos((angle + M_PI) / 2)), 2);
 	}
 
 	for (size_t i=0; i<=N; i++) {
diff --git a/hspm/src/solver.cpp b/hspm/src/solver.cpp
index 4e986a9..8972040 100644
--- a/hspm/src/solver.cpp
+++ b/hspm/src/solver.cpp
@@ -34,7 +34,7 @@ void HSPM::solve()
     // q = s1*tau + s2
     // We need Kutta for tau
     tau = this->solveOffBodyKutta();
-    
+
     q = s1*tau + s2;
 
     this->computeInviscidVelocity(); 
@@ -84,6 +84,9 @@ double HSPM::solveOffBodyKutta() {
     double _tau_terms = 2*(_a*_b + _c*_d - _e*_f - _g*_h);
     double _const_terms = pow(_b,2) + pow(_d,2) - pow(_f,2) - pow(_h,2);
 
+    if (_tau2_terms < 1e-10) 
+        return -_const_terms / _tau_terms; // If the quadratic equation is degenerate, we can solve it directly
+
     // Solve the quadratic equation
     double tau_ = (-_tau_terms + sqrt(pow(_tau_terms,2) - 4*_tau2_terms*_const_terms)) / (2*_tau2_terms);
 
diff --git a/python/main.py b/python/main.py
index 44e3b17..d8d9cd1 100644
--- a/python/main.py
+++ b/python/main.py
@@ -1,23 +1,16 @@
 from utils import *
 
-N = 1000
-
 config = {
 	'chord': 1,
-	'it_solver_tolerance': 1e-6,
 	'aoa': 2,
 	'naca': "0012",
-	'N': N,
+	'N': 100,
 	'filename': "airfoil.dat"
 }
 
 if __name__ == "__main__":
 	solver = initHSPM(config)
 
-	"""for i in range(N):
-					solver.imposeBlowingVelocity(i, 1)
-					solver.setdStar(i, 0.1)"""
-
 	solver.solve()
 
 	cl = solver.cl
-- 
GitLab