cleanup of massActions...and change definition of K in the momas example

Author: rrsettgast

src/reactions/geochemistry/Momas.hpp

@@ -38,13 +38,13 @@ constexpr CArrayWrapper<double, 7, 12> stoichMatrix =
 
 constexpr CArrayWrapper<double, 7> equilibriumConstants = 
   { 
-    1.0e-12,   // C1 = -X2
-    1.0,       // C2 = X2 + X3
-    1.0,       // C3 = X2 + X4
-    0.1,       // C4 = -4X2 + X3 + 3X4
-    1.0e35,    // C5 = 4X2 + 3X3 + X4
-    1.0e6,     // CS1 = 3X2 + X3 + S
-    1.0e-1     // CS2 = -3X2 + X4 + 2S
+    1.0e12,   //   C1 + X2 = ??
+    1.0,      //        C2 = X2 + X3
+    1.0,      //        C3 = X2 + X4
+    1.0e1,    //  C4 + 4X2 = X3 + 3X4
+    1.0e-35,  //        C5 = 4X2 + 3X3 + X4
+    1.0e-6,   //       CS1 = 3X2 + X3 + S
+    1.0e1     // CS2 + 3X2 = + X4 + 2S
   };
 
 constexpr CArrayWrapper<double, 7> forwardRates = 

src/reactions/massActions/MassActions.hpp

@@ -172,42 +172,17 @@ void calculateAggregatePrimaryConcentrationsWrtLogC( PARAMS_DATA const & params,
                                                      ARRAY_2D & dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations )
 {
   constexpr int numSecondarySpecies = PARAMS_DATA::numSecondarySpecies();
-  constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();
 
   REAL_TYPE logSecondarySpeciesConcentrations[numSecondarySpecies] = {0};
 
-  for( INDEX_TYPE i = 0; i < numPrimarySpecies; ++i )
-  {
-    for( INDEX_TYPE j = 0; j < numPrimarySpecies; ++j )
-    {
-      dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] = 0.0;
-    }
-  }
-
-  calculateLogSecondarySpeciesConcentration< REAL_TYPE,
-                                             INT_TYPE,
-                                             INDEX_TYPE >( params,
-                                                           logPrimarySpeciesConcentrations,
-                                                           logSecondarySpeciesConcentrations );
+  calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, 
+                                                  INT_TYPE, 
+                                                  INDEX_TYPE>( params,
+                                                               logPrimarySpeciesConcentrations,
+                                                               logSecondarySpeciesConcentrations,
+                                                               aggregatePrimarySpeciesConcentrations,
+                                                               dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations );
 
-  for( int i = 0; i < numPrimarySpecies; ++i )
-  {
-    REAL_TYPE const speciesConcentration_i = exp( logPrimarySpeciesConcentrations[i] );
-    aggregatePrimarySpeciesConcentrations[i] = speciesConcentration_i;
-    dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, i ) = speciesConcentration_i;
-    for( int j = 0; j < numSecondarySpecies; ++j )
-    {
-      REAL_TYPE const secondarySpeciesConcentrations_j = exp( logSecondarySpeciesConcentrations[j] );
-      aggregatePrimarySpeciesConcentrations[i] += params.stoichiometricMatrix( j, i+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
-      for( int k=0; k<numPrimarySpecies; ++k )
-      {
-        REAL_TYPE const dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration = params.stoichiometricMatrix( j, k+numSecondarySpecies ) * secondarySpeciesConcentrations_j;
-        dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations( i, k ) += params.stoichiometricMatrix( j,
-                                                                                                                                   i+numSecondarySpecies ) *
-                                                                                                      dSecondarySpeciesConcentrations_dLogPrimarySpeciesConcentration;
-      }
-    }
-  }
 }
 
 template< typename REAL_TYPE,

src/reactions/reactionsSystems/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp

@@ -95,6 +95,9 @@ EquilibriumReactions< REAL_TYPE,
 
 
   REAL_TYPE residualNorm = 0.0;
+ //         0:     1e-20       -0           2 -2.5e+11       1e-20        7           2      1.8           1        5  
+  printf( "iter       X1       R0           X2      R1          X3       R2          X4       R3           S       R4\n");
+  printf( "----   ---------------      ---------------      ---------------      ---------------      ---------------\n");
   for( int k=0; k<30; ++k )
   {
     computeResidualAndJacobianAggregatePrimaryConcentrations( temperature,
@@ -110,7 +113,21 @@ EquilibriumReactions< REAL_TYPE,
       residualNorm += residual[i] * residual[i];
     }
     residualNorm = sqrt( residualNorm );
-    printf( "iter, residualNorm = %2d, %16.10g \n", k, residualNorm );
+
+      printf( "%2d:  %8.2g %8.2g    %8.2g %8.2g    %8.2g %8.2g    %8.2g %8.2g    %8.2g %8.2g  \n", 
+              k,
+              exp( logPrimarySpeciesConcentration[0] ),
+              residual[0],
+              exp( logPrimarySpeciesConcentration[1] ),
+              residual[1],
+              exp( logPrimarySpeciesConcentration[2] ),
+              residual[2],
+              exp( logPrimarySpeciesConcentration[3] ),
+              residual[3],
+              exp( logPrimarySpeciesConcentration[4] ),
+              residual[4] );
+
+    //printf( "iter, residualNorm = %2d, %16.10g \n", k, residualNorm );
     if( residualNorm < 1.0e-12 )
     {
       printf( " converged\n" );