more testing of kinetic reactions. using ln(c) results in a very poorly conditioned system

Author: rrsettgast

src/common/DirectSystemSolve.hpp

@@ -7,55 +7,59 @@ namespace hpcReact
 template< typename REAL_TYPE, int N >
 HPCREACT_HOST_DEVICE void solveNxN_pivoted(REAL_TYPE A[N][N], REAL_TYPE b[N], REAL_TYPE x[N]) 
 {
-    int pivot[N] = {0, 1, 2};  // Row index tracker
-    for (int i = 0; i < N; i++) 
+
+  
+  int pivot[N] = {0, 1, 2};  // Row index tracker
+  for (int i = 0; i < N; i++) 
+  {
+    pivot[i] = i;
+  }
+
+  // **Step 1: Forward Elimination with Pivoting**
+  for (int k = 0; k < N-1; k++) 
+  {
+    // **Find Pivot Row**
+    int max_row = k;
+    REAL_TYPE max_val = abs(A[pivot[k]][k]);
+    for (int i = k + 1; i < N; i++) 
     {
-        pivot[i] = i;
+      if (fabs(A[pivot[i]][k]) > max_val) 
+      {
+        max_val = fabs(A[pivot[i]][k]);
+        max_row = i;
+      }
     }
 
-    // **Step 1: Forward Elimination with Pivoting**
-    for (int k = 0; k < N-1; k++) 
+    // **Swap Rows in Pivot Array**
+    if (max_row != k) 
     {
-        // **Find Pivot Row**
-        int max_row = k;
-        REAL_TYPE max_val = abs(A[pivot[k]][k]);
-        for (int i = k + 1; i < N; i++) 
-        {
-            if (fabs(A[pivot[i]][k]) > max_val) 
-            {
-                max_val = fabs(A[pivot[i]][k]);
-                max_row = i;
-            }
-        }
-
-        // **Swap Rows in Pivot Array**
-        if (max_row != k) 
-        {
-            int temp = pivot[k];
-            pivot[k] = pivot[max_row];
-            pivot[max_row] = temp;
-        }
-
-        // **Gaussian Elimination**
-        for (int i = k + 1; i < N; i++) 
-        {
-            REAL_TYPE factor = A[pivot[i]][k] / A[pivot[k]][k];
-            for (int j = k; j < N; j++) 
-            {
-                A[pivot[i]][j] -= factor * A[pivot[k]][j];
-            }
-            b[pivot[i]] -= factor * b[pivot[k]];
-        }
+      int temp = pivot[k];
+      pivot[k] = pivot[max_row];
+      pivot[max_row] = temp;
     }
 
-    // **Step 2: Back-Substitution**
-    for (int i = N - 1; i >= 0; --i) {
-        x[i] = b[pivot[i]];
-        for (int j = i + 1; j < N; j++) {
-            x[i] -= A[pivot[i]][j] * x[j];
-        }
-        x[i] /= A[pivot[i]][i]; // Normalize
-}
+    // **Gaussian Elimination**
+    for (int i = k + 1; i < N; i++) 
+    {
+      REAL_TYPE factor = A[pivot[i]][k] / A[pivot[k]][k];
+      for (int j = k; j < N; j++) 
+      {
+        A[pivot[i]][j] -= factor * A[pivot[k]][j];
+      }
+      b[pivot[i]] -= factor * b[pivot[k]];
+    }
+  }
+
+  // **Step 2: Back-Substitution**
+  for (int i = N - 1; i >= 0; --i) 
+  {
+    x[i] = b[pivot[i]];
+    for (int j = i + 1; j < N; j++) 
+    {
+      x[i] -= A[pivot[i]][j] * x[j];
+    }
+    x[i] /= A[pivot[i]][i]; // Normalize
+  }
 
 }
 

src/reactions/bulkGeneric/KineticReactions_impl.hpp

@@ -256,15 +256,13 @@ HPCREACT_HOST_DEVICE inline void
 KineticReactions< REAL_TYPE,
                   INT_TYPE,
                   INDEX_TYPE,
-                  LOGE_CONCENTRATION
-                  >::timeStep( 
-                    RealType const dt,
-                    RealType const & temperature,
-                    PARAMS_DATA const & params,
-                    ARRAY_1D_TO_CONST const & speciesConcentration_n,
-                    ARRAY_1D & speciesConcentration,
-                    ARRAY_1D & speciesRates,
-                    ARRAY_2D & speciesRatesDerivatives )
+                  LOGE_CONCENTRATION >::timeStep( RealType const dt,
+                                                  RealType const & temperature,
+                                                  PARAMS_DATA const & params,
+                                                  ARRAY_1D_TO_CONST const & speciesConcentration_n,
+                                                  ARRAY_1D & speciesConcentration,
+                                                  ARRAY_1D & speciesRates,
+                                                  ARRAY_2D & speciesRatesDerivatives )
 {
 //  constexpr int numReactions = PARAMS_DATA::numReactions;
   constexpr int numSpecies = PARAMS_DATA::numSpecies;
@@ -285,32 +283,82 @@ KineticReactions< REAL_TYPE,
 
     double residual[numSpecies] = { 0.0 };
     double deltaPrimarySpeciesConcentration[numSpecies] = { 0.0 };
+
+    // form residual and Jacobian
     for( int i = 0; i < numSpecies; ++i )
     {
-      residual[i] = -(speciesConcentration[i] - speciesConcentration_n[i] - dt * speciesRates[i]);
+
+
+      RealType nonLogC;
+      RealType nonLogC_n;
+      if constexpr( LOGE_CONCENTRATION )
+      {
+        nonLogC = exp( speciesConcentration[i] );
+        nonLogC_n = exp( speciesConcentration_n[i] );
+      }
+      else
+      {
+        nonLogC = speciesConcentration[i];
+        nonLogC_n = speciesConcentration_n[i];
+      }
+      residual[i] = -(nonLogC - nonLogC_n - dt * speciesRates[i]);
+
+
       for( int j = 0; j < numSpecies; ++j )
       {
         speciesRatesDerivatives( i, j ) = - dt * speciesRatesDerivatives( i, j );
       }
-      speciesRatesDerivatives( i, i ) += 1.0;
+      if constexpr( LOGE_CONCENTRATION )
+      {
+        speciesRatesDerivatives( i, i ) += nonLogC;
+      }
+      else
+      {
+        speciesRatesDerivatives( i, i ) += 1.0;
+      }
     }
 
+
+
     residualNorm = 0.0;
     for( int j = 0; j < numSpecies; ++j )
     {
       residualNorm += residual[j] * residual[j];
     }
     residualNorm = sqrt( residualNorm );
-    if( residualNorm < 1.0e-8 )
+    if( residualNorm < 1.0e-14 )
     {
       break;
     }
 
+
+//     printf( "residual = { ");
+//     for( int i = 0; i < numSpecies; ++i )
+//     {
+//       printf( " %g, ", residual[i] );
+//     }
+//     printf( "}\n" );
+
+//     printf( "Jacobian = { \n" );
+//     for( int i = 0; i < numSpecies; ++i )
+//     {
+//       printf( " { ");
+//       for( int j = 0; j < numSpecies; ++j )
+//       {
+//         printf( " %g ", speciesRatesDerivatives( i, j ) );
+// //        printf( " %g ", speciesRatesDerivatives( i, j ) / exp(speciesConcentration[j]) );
+//         if( j < numSpecies-1 ) { printf( ", " ); }
+//       }
+//       printf( "}, \n" );
+//     }
+//     printf( "}\n" );
+
     solveNxN_pivoted<double,numSpecies>( speciesRatesDerivatives.data, residual, deltaPrimarySpeciesConcentration );
 
     for( int i = 0; i < numSpecies; ++i )
     {
-      speciesConcentration[i] += deltaPrimarySpeciesConcentration[i];
+//      printf( "species %2d: concentration = %e, residual = %e, delta = %e \n", i, speciesConcentration[i], residual[i], deltaPrimarySpeciesConcentration[i] );
+      speciesConcentration[i] = speciesConcentration[i] + deltaPrimarySpeciesConcentration[i];
     }
 
   }