fixed more floating point equality comparisons & adjusted tolerance f…

…or floating point equality in 1d point control
Cantera · Jun 23, 2023 · 00963a6 · 00963a6
1 parent 3b195b8
commit 00963a6
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 10 deletions.
diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp
@@ -753,8 +753,7 @@ double Sim1D::fixedTemperatureLocation()
 
 void Sim1D::setFuelInternalBoundary(double temperature)
 {
-    double epsilon = 1e-10; // Define your precision threshold here
-
+    double epsilon = 1e-5; // Precision threshold for being 'equal' to a temperature
     for (size_t n = 0; n < nDomains(); n++) {
         Domain1D& d = domain(n);
 

diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp
@@ -404,7 +404,6 @@ void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax)
 
 void StFlow::computeRadiation(double* x, size_t jmin, size_t jmax) 
 {
-
     // calculation of qdotRadiation
 
     // The simple radiation model used was established by Y. Liu and B. Rogg [Y.
@@ -579,8 +578,8 @@ void StFlow::evalRadialMomentumEquation(double* x, double* rsd, int* diag, doubl
 
 void StFlow::evalContinuityEquation(double* x, double* rsd, int* diag, double rdt, size_t j) 
 {
-    //algebraic constraint
-    diag[index(c_offset_U, j)] = 0;
+
+    double epsilon = 1e-5; // Precision threshold for being 'equal' to a coordinate
     //----------------------------------------------
     //    Continuity equation
     //
@@ -598,7 +597,7 @@ void StFlow::evalContinuityEquation(double* x, double* rsd, int* diag, double rd
             rsd[index(c_offset_U,j)] =
                 - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]
                 - (density(j-1)*V(x,j-1) + density(j)*V(x,j));
-        } else if (grid(j) == m_zfixed) {
+        } else if (std::abs(grid(j) - m_zfixed) < epsilon) {
             if (m_do_energy[j]) {
                 rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
             } else {
@@ -611,6 +610,8 @@ void StFlow::evalContinuityEquation(double* x, double* rsd, int* diag, double rd
                 - (density(j+1)*V(x,j+1) + density(j)*V(x,j));
         }
     }
+    //algebraic constraint
+    diag[index(c_offset_U, j)] = 0;
 }
 
 void StFlow::evalSpeciesEquation(double* x, double* rsd, int* diag, double rdt, size_t j)
@@ -676,18 +677,19 @@ void StFlow::evalEnergyEquation(double* x, double* rsd, int* diag, double rdt, s
 
 void StFlow::evalLEquation(double* x, double* rsd, int* diag, double rdt, size_t j) 
 {
+    double epsilon = 1e-5; // Precision threshold for being 'equal' to a coordinate
     if (m_onePointControl) {
         if (grid(j) > m_zFuel) {
             rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
-        } else if (grid(j) == m_zFuel) {
+        } else if (std::abs(grid(j) - m_zFuel) < epsilon ) {
             rsd[index(c_offset_L, j)] = T(x,j) - m_tFuel;
         } else if (grid(j) < m_zFuel) {
             rsd[index(c_offset_L, j)] = lambda(x,j+1) - lambda(x,j);
         }
     } else if (m_twoPointControl) {
         if (grid(j) > m_zFuel) {
             rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
-        } else if (grid(j) == m_zFuel) {
+        } else if (std::abs(grid(j) - m_zFuel) < epsilon ) {
             rsd[index(c_offset_L, j)] = T(x,j) - m_tFuel;
         } else if (grid(j) < m_zFuel) {
             rsd[index(c_offset_L, j)] = lambda(x,j+1) - lambda(x,j);
@@ -700,12 +702,13 @@ void StFlow::evalLEquation(double* x, double* rsd, int* diag, double rdt, size_t
 
 void StFlow::evalUoEquation(double* x, double* rsd, int* diag, double rdt, size_t j) 
 {
+    double epsilon = 1e-5; // Precision threshold for being 'equal' to a coordinate
     if (m_onePointControl) {
         rsd[index(c_offset_Uo, j)] = Uo(x,j) - Uo(x,j-1);
     } else if (m_twoPointControl) {
         if (grid(j) > m_zOxid) {
             rsd[index(c_offset_Uo, j)] = Uo(x,j) - Uo(x,j-1);
-        } else if (grid(j) == m_zOxid) {
+        } else if (std::abs(grid(j) - m_zOxid) < epsilon) {
             rsd[index(c_offset_Uo, j)] = T(x,j) - m_tOxid;
         } else if (grid(j) < m_zOxid) {
             rsd[index(c_offset_Uo, j)] = Uo(x,j+1) - Uo(x,j);
@@ -1105,7 +1108,7 @@ void StFlow::setMeta(const AnyMap& state)
             m_twoPointControl = false;
             m_zFuel = pc["location"].asDouble();
             m_tFuel = pc["temperature"].asDouble();
-        } else if (state["fixed-point"]["type"] == "two-point") {
+        } else if (pc["type"] == "two-point") {
             m_onePointControl = false;
             m_twoPointControl = true;
             m_zFuel = pc["fuel-location"].asDouble();