eps change

bigfooted · bigfooted · commit 04f0681517b3 · 2025-12-21T20:14:14.000+01:00
diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp
@@ -505,15 +505,17 @@ void CSysMatrix<ScalarType>::Gauss_Elimination(ScalarType* matrix, ScalarType* v
 #define A(I, J) matrix[(I)*nVar + (J)]
 
   /*--- Regularization epsilon to prevent divide-by-zero ---*/
-  constexpr passivedouble eps = 1e-12;
+  //constexpr passivedouble eps = 1e-12;
+  const su2double eps = 1e-12;
 
   /*--- Transform system in Upper Matrix ---*/
 
   for (auto iVar = 1ul; iVar < nVar; iVar++) {
     for (auto jVar = 0ul; jVar < iVar; jVar++) {
       /*--- Regularize pivot if too small to prevent divide-by-zero ---*/
       if (std::abs(A(jVar, jVar)) < eps) {
-        A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps);
+        //A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps);
+        A(jVar, jVar) = (A(jVar, jVar) >= 0.0 ? eps : -eps);
         std::cout << "DEBUG Gauss_Elimination: Regularized small pivot A(" << jVar << "," << jVar << ") to "
                   << A(jVar, jVar) << std::endl;
       }
@@ -532,7 +534,8 @@ void CSysMatrix<ScalarType>::Gauss_Elimination(ScalarType* matrix, ScalarType* v
 
     /*--- Regularize diagonal if too small ---*/
     if (std::abs(A(iVar, iVar)) < eps) {
-      A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps);
+      //A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps);
+      A(iVar, iVar) = (A(iVar, iVar) >= 0.0 ? eps : -eps);
       std::cout << "DEBUG Gauss_Elimination backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to "
                 << A(iVar, iVar) << std::endl;
     }
@@ -568,14 +571,14 @@ void CSysMatrix<ScalarType>::MatrixInverse(ScalarType* matrix, ScalarType* inver
 #define A(I, J) matrix[(I)*nVar + (J)]
 
   /*--- Regularization epsilon to prevent divide-by-zero ---*/
-  constexpr passivedouble eps = 1e-12;
+  const su2double eps = 1e-12;
 
   /*--- Transform system in Upper Matrix ---*/
   for (auto iVar = 1ul; iVar < nVar; iVar++) {
     for (auto jVar = 0ul; jVar < iVar; jVar++) {
       /*--- Regularize pivot if too small to prevent divide-by-zero ---*/
       if (std::abs(A(jVar, jVar)) < eps) {
-        A(jVar, jVar) = (A(jVar, jVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps);
+        A(jVar, jVar) = (A(jVar, jVar) >= 0.0 ? eps : -eps);
         std::cout << "DEBUG MatrixInverse: Regularized small pivot A(" << jVar << "," << jVar << ") to "
                   << A(jVar, jVar) << " at iVar=" << iVar << std::endl;
       }
@@ -596,7 +599,7 @@ void CSysMatrix<ScalarType>::MatrixInverse(ScalarType* matrix, ScalarType* inver
 
     /*--- Regularize diagonal if too small ---*/
     if (std::abs(A(iVar, iVar)) < eps) {
-      A(iVar, iVar) = (A(iVar, iVar) >= ScalarType(0)) ? ScalarType(eps) : ScalarType(-eps);
+      A(iVar, iVar) = (A(iVar, iVar) >= 0.0 ? eps : -eps);
       std::cout << "DEBUG MatrixInverse backsubst: Regularized small diagonal A(" << iVar << "," << iVar << ") to "
                 << A(iVar, iVar) << std::endl;
     }
diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/QUICK_FIX_GUIDE.md b/TestCases/py_wrapper/turbulent_premixed_psi/QUICK_FIX_GUIDE.md
@@ -0,0 +1,124 @@
+# Quick Fix Summary: Residual Stagnation Solution
+
+## The Problem
+```
+Species equation → converges ✓
+Update T from species → T = (1-c)*Tu + c*Tf
+Flow equations → STAGNATE ✗
+```
+
+**Why?** Density wasn't updated when temperature changed!
+
+## The Fix
+
+### BEFORE (Broken):
+```python
+def update_temperature(SU2Driver, iPoint):
+    C = SU2Driver.Solution(iSPECIESSOLVER)(iPoint,0)
+    T = Tu*(1-C) + Tf*C
+    Cp = SU2Driver.Primitives()(iPoint,iCp)
+    SU2Driver.Solution(iFLOWSOLVER).Set(iPoint,iTEMP,Cp*T)
+    # ⚠️ Density not updated! ρ still has old value
+```
+
+### AFTER (Fixed):
+```python
+def update_temperature(SU2Driver, iPoint):
+    C = SU2Driver.Solution(iSPECIESSOLVER)(iPoint,0)
+    T_target = Tu*(1-C) + Tf*C
+    T_old = SU2Driver.Primitives()(iPoint,iTEMPERATURE)
+
+    # ✓ Under-relaxation for stability
+    T_new = RELAX_FACTOR * T_target + (1.0 - RELAX_FACTOR) * T_old
+
+    # ✓ Compute gas constant from old state
+    R_gas = P / (rho_old * T_old)
+
+    # ✓ Update density consistently: ρ = P/(RT)
+    rho_new = P / (R_gas * T_new)
+
+    # ✓ Update both temperature and density
+    SU2Driver.Solution(iFLOWSOLVER).Set(iPoint, iTEMP, T_new)
+    SU2Driver.Primitives().Set(iPoint, iDENSITY, rho_new)
+```
+
+## Main Loop Sequence
+
+### BEFORE:
+```python
+driver.Preprocess()
+driver.Run()
+# Set source terms (WRONG: after Run!)
+for i in nodes: Source.Set(i, zimont(i))
+# Update temperature (no density update)
+for i in nodes: update_temperature(i)
+driver.Postprocess()
+```
+
+### AFTER:
+```python
+driver.Preprocess()
+# ✓ Set source terms BEFORE Run
+for i in nodes: Source.Set(i, zimont(i))
+# ✓ Run solvers
+driver.Run()
+# ✓ Update temperature AND density AFTER Run
+for i in nodes: update_temperature(i)
+# ✓ Monitor coupling
+print(f"Max ΔT: {max_delta_T}, Max Δρ/ρ: {max_delta_rho}")
+driver.Postprocess()
+```
+
+## Key Parameters
+
+```python
+# Add at top of file:
+RELAX_FACTOR = 0.7  # Adjust 0.5-0.8 for stability vs speed
+```
+
+## Expected Behavior
+
+| Before Fix | After Fix |
+|------------|-----------|
+| Species: ✓ converging | Species: ✓ converging |
+| Pressure: ✗ stagnating | Pressure: ✓ converging |
+| Velocity: ✗ stagnating | Velocity: ✓ converging |
+| ρ inconsistent with T | ρ = P/(RT) ✓ |
+
+## Physics
+
+**In combustion:**
+- c: 0 → 1 (progress)
+- T: 673K → 1777K ↑
+- ρ: 2.52 → ~0.95 ↓ (MUST decrease!)
+
+**If ρ not updated:**
+- Continuity: wrong mass fluxes
+- Momentum: wrong inertia
+- → Residuals can't converge!
+
+## Monitoring
+
+Look for in output:
+```
+Source term - Max: 1.23e+02, Min: 0.00e+00, Avg: 3.45e+01
+Max temperature change: 5.23e+01 K
+Max relative density change: 2.15e-02
+```
+
+## Troubleshooting
+
+**Still not converging?**
+1. Reduce `RELAX_FACTOR` to 0.5
+2. Reduce CFL number in config
+3. Check source term isn't too strong
+4. Verify mass conservation: Σ(source terms) ≈ 0
+
+**Oscillations?**
+1. Reduce `RELAX_FACTOR`
+2. Check grid quality in flame region
+3. Ensure boundary conditions are correct
+
+## One-Line Summary
+
+**Always update density when temperature changes: `ρ_new = P/(R*T_new)`**
diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/SOLUTION_NOTES.md b/TestCases/py_wrapper/turbulent_premixed_psi/SOLUTION_NOTES.md
@@ -0,0 +1,128 @@
+# Solution for Residual Stagnation in Turbulent Premixed Combustion
+
+## Problem Description
+
+When using the Python wrapper to add species transport with strong source terms and algebraically updating temperature/enthalpy from species composition, the pressure and velocity residuals stagnate in regions with high source terms while the species residuals converge properly.
+
+## Root Cause
+
+The issue occurs because:
+
+1. **Species solver** converges and updates the progress variable `c`
+2. **Temperature is updated algebraically**: `T = (1-c)*Tu + c*Tf`
+3. **Density is NOT updated** to match the new temperature
+4. **Thermodynamic inconsistency**: The ideal gas law `ρ = P/(RT)` is violated
+5. **Flow equations see inconsistent state**: Continuity and momentum equations operate with outdated density
+
+This breaks the coupling between species, energy, and momentum equations, causing residual stagnation.
+
+## Solution Implemented (Solution #1)
+
+### Key Changes in `run.py`:
+
+#### 1. **Updated `update_temperature()` function** (Lines ~90-130)
+   - Now computes new density when temperature changes: `ρ_new = P/(R*T_new)`
+   - Updates both temperature AND density in the primitive variables
+   - Added under-relaxation factor for stability: `T_new = ω*T_target + (1-ω)*T_old`
+   - Maintains thermodynamic consistency with ideal gas law
+
+#### 2. **Improved iteration sequence** (Lines ~250-290)
+   - **STEP 1**: Set species source terms
+   - **STEP 2**: Run all solvers (species equation solved)
+   - **STEP 3**: Update temperature AND density from new species field
+   - **STEP 4**: Postprocess and update
+   - Ensures proper coupling between species and flow equations
+
+#### 3. **Added under-relaxation** (Line ~58)
+   - `RELAX_FACTOR = 0.7` prevents oscillations
+   - Can be tuned between 0.5-0.8 depending on source term strength
+   - Lower values = more stable but slower convergence
+
+#### 4. **Added coupling diagnostics** (Lines ~270-280)
+   - Monitors `max_delta_T` and `max_delta_rho`
+   - Helps identify coupling issues
+   - Shows when thermodynamic updates are stabilizing
+
+#### 5. **Added source term monitoring** (Lines ~185-215, ~260-265)
+   - Tracks max/min/avg source term values
+   - Helps diagnose convergence problems
+   - Important for checking mass conservation
+
+## Physical Interpretation
+
+For a premixed flame:
+- Progress variable goes from `c=0` (unburnt) to `c=1` (burnt)
+- Temperature increases: `T: 673K → 1777K`
+- Density MUST decrease (isobaric): `ρ ∝ 1/T`
+- Without density update: continuity equation sees wrong mass fluxes
+- Result: momentum residuals can't converge
+
+## Expected Results
+
+After this fix:
+- ✅ Species residuals converge (as before)
+- ✅ Temperature correctly updates from species
+- ✅ **Density correctly updates from temperature**
+- ✅ **Pressure residuals converge** (no longer stagnate)
+- ✅ **Velocity residuals converge** (no longer stagnate)
+- ✅ Thermodynamic consistency maintained
+
+## Tuning Parameters
+
+If convergence is still slow, adjust:
+
+### 1. **Under-relaxation factor** (Line 58)
+```python
+RELAX_FACTOR = 0.5-0.8  # Lower = more stable, slower
+```
+
+### 2. **CFL number in config file**
+```
+CFL_NUMBER= 10.0  # Reduce if needed
+CFL_ADAPT= YES
+CFL_ADAPT_PARAM= ( 0.3, 0.5, 10.0, 100.0, 1e-6 )
+```
+
+### 3. **Number of inner iterations**
+```python
+for inner_iter in range(10):  # Increase for more coupling iterations
+```
+
+## Additional Recommendations
+
+1. **Monitor residual ratios**: Watch for regions where flow residuals are high relative to species residuals
+
+2. **Check mass conservation**: The source term should not create/destroy mass:
+   ```
+   ∂(ρc)/∂t + ∇·(ρcv) = S_c
+   ```
+   where `S_c` is per unit volume
+
+3. **Verify ideal gas law**: At any point, check if `P = ρRT` holds after updates
+
+4. **Use implicit time stepping**: Helps with stiff source terms (already in config)
+
+5. **Consider local time stepping**: Can help in regions with varying source terms
+
+## Technical Details
+
+### Incompressible Flow Considerations
+This case uses `INC.FLOW` (incompressible solver), which:
+- Solves for pressure directly (not density-based)
+- Assumes low Mach number
+- Still requires density updates for variable density flows (like combustion!)
+- Energy equation is decoupled but affects momentum through density
+
+### Variable Density Incompressible Flow
+Even though it's "incompressible", combustion creates variable density:
+- ∇·(ρv) = 0 (mass conservation, not ∇·v = 0)
+- ρ varies due to temperature, not pressure waves
+- Density MUST be updated when temperature changes
+
+## References
+- See main SU2 documentation for species transport
+- Zimont turbulent flame speed closure: Zimont, V.L., 2000. "Gas premixed combustion at high turbulence."
+- Variable density incompressible flow formulation
+
+## Contact
+For issues or questions about this fix, refer to the SU2 forums or GitHub issues.
diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/VISUALIZATION.md b/TestCases/py_wrapper/turbulent_premixed_psi/VISUALIZATION.md
