Skip to content

Commit 50587d6

Browse files
rrsettgastrrsettgast
committed
fixed bug
1 parent 504d07a commit 50587d6

File tree

4 files changed

+91
-59
lines changed

4 files changed

+91
-59
lines changed

‎src/common/DirectSystemSolve.hpp‎

Lines changed: 14 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ bool isPositiveDefinite( REAL_TYPE const (&A)[N][N] )
2525
{
2626
// Compute determinant of k-th leading principal minor using Gaussian elimination
2727
if( temp[k][k] <= 0 )
28-
return false; // Must be positive
28+
return false; // Must be positive
2929

3030
for( int i = k + 1; i < N; i++ )
3131
{
@@ -63,19 +63,19 @@ void solveNxN_Cholesky( REAL_TYPE const (&A)[N][N], REAL_TYPE const (&b)[N], REA
6363
}
6464

6565
// **Forward Substitution: Solve L y = b**
66-
REAL_TYPE y[N];
66+
//REAL_TYPE y[N];
6767
for( int i = 0; i < N; i++ )
6868
{
69-
y[i] = b[i];
69+
x[i] = b[i];
7070
for( int j = 0; j < i; j++ )
71-
y[i] -= L[i][j] * y[j];
72-
y[i] /= L[i][i];
71+
x[i] -= L[i][j] * x[j];
72+
x[i] /= L[i][i];
7373
}
7474

7575
// **Backward Substitution: Solve L^T x = y**
7676
for( int i = N - 1; i >= 0; i-- )
7777
{
78-
x[i] = y[i];
78+
// x[i] = y[i];
7979
for( int j = i + 1; j < N; j++ )
8080
x[i] -= L[j][i] * x[j];
8181
x[i] /= L[i][i];
@@ -102,9 +102,13 @@ void solveNxN_Cholesky( symmetricMatrix< REAL_TYPE, int, N > const & A,
102102
sum += L( i, k ) * L( j, k );
103103
}
104104
if( i == j )
105+
{
105106
L( i, j ) = sqrt( A( i, i ) - sum );
107+
}
106108
else
109+
{
107110
L( i, j ) = (A( i, j ) - sum) / L[j][j];
111+
}
108112
}
109113
}
110114

@@ -114,7 +118,9 @@ void solveNxN_Cholesky( symmetricMatrix< REAL_TYPE, int, N > const & A,
114118
{
115119
y[i] = b[i];
116120
for( int j = 0; j < i; j++ )
121+
{
117122
y[i] -= L( i, j ) * y[j];
123+
}
118124
y[i] /= L( i, i );
119125
}
120126

@@ -123,7 +129,9 @@ void solveNxN_Cholesky( symmetricMatrix< REAL_TYPE, int, N > const & A,
123129
{
124130
x[i] = y[i];
125131
for( int j = i + 1; j < N; j++ )
132+
{
126133
x[i] -= L( j, i ) * x[j];
134+
}
127135
x[i] /= L( i, i );
128136
}
129137
}

‎src/common/symmetricMatrix.hpp‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
#pragma once
2+
3+
template< typename T, typename INDEX_TYPE, int N >
4+
struct symmetricMatrix
5+
{
6+
static constexpr INDEX_TYPE size() { return ( N*(N+1) ) / 2; }
7+
8+
static inline INDEX_TYPE linearIndex( INDEX_TYPE const i, INDEX_TYPE const j )
9+
{
10+
return (( i*(i+1) ) >> 1) + j;
11+
}
12+
13+
inline T & operator()( INDEX_TYPE const i, INDEX_TYPE const j )
14+
{
15+
return m_data[linearIndex(i,j)];
16+
}
17+
18+
inline T const & operator()( INDEX_TYPE const i, INDEX_TYPE const j ) const
19+
{
20+
return m_data[linearIndex(i,j)];
21+
}
22+
23+
T m_data[ size() ];
24+
25+
};

‎src/reactions/bulkGeneric/EquilibriumReactions_impl.hpp‎

Lines changed: 33 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,9 @@ EquilibriumReactions< REAL_TYPE,
6868

6969
RealType forwardProduct = 1.0;
7070
RealType reverseProduct = 1.0;
71+
72+
// these actually only hold the derivatives of the single concentration term for the product...not the full product.
73+
// it will have to be multiplied by the product itself to get the derivative of the product.
7174
RealType dForwardProduct_dxi[numReactions] = {0.0};
7275
RealType dReverseProduct_dxi[numReactions] = {0.0};
7376
// loop over species
@@ -98,15 +101,15 @@ EquilibriumReactions< REAL_TYPE,
98101
// compute the residual for this reaction
99102
if constexpr( RESIDUAL_FORM == 0 )
100103
{
101-
residual[a] = forwardProduct - Keq * reverseProduct;
104+
residual[a] = Keq * forwardProduct - reverseProduct;
102105
}
103106
else if constexpr( RESIDUAL_FORM == 1 )
104107
{
105-
residual[a] = 1.0 - Keq * reverseProduct / forwardProduct;
108+
residual[a] = 1.0 - reverseProduct / ( forwardProduct * Keq );
106109
}
107110
else if constexpr( RESIDUAL_FORM == 2 )
108111
{
109-
residual[a] = log( Keq * reverseProduct / forwardProduct );
112+
residual[a] = log( reverseProduct / ( forwardProduct * Keq ) );
110113
}
111114
// printf( "residual[%d] = %g - %g * %g = %g\n", a, forwardProduct, Keq, reverseProduct, residual[a] );
112115

@@ -117,20 +120,25 @@ EquilibriumReactions< REAL_TYPE,
117120
// printf( "(%d) dReverseProduct_dxi[%d] = %f\n", a, b, dReverseProduct_dxi[b] );
118121
// printf( "Keq = %f\n", Keq );
119122

120-
dForwardProduct_dxi[b] *= forwardProduct;
121-
dReverseProduct_dxi[b] *= reverseProduct;
122123

123124
if constexpr( RESIDUAL_FORM == 0 )
124125
{
125-
jacobian( a, b ) = dForwardProduct_dxi[b] - Keq * dReverseProduct_dxi[b];
126+
dForwardProduct_dxi[b] *= forwardProduct;
127+
dReverseProduct_dxi[b] *= reverseProduct;
128+
jacobian( a, b ) = Keq * dForwardProduct_dxi[b] - dReverseProduct_dxi[b];
126129
}
127130
else if constexpr( RESIDUAL_FORM == 1 )
128131
{
129-
jacobian( a, b ) = -Keq * ( dReverseProduct_dxi[b] / forwardProduct - dForwardProduct_dxi[b] * reverseProduct / ( forwardProduct * forwardProduct ) );
132+
dForwardProduct_dxi[b] *= forwardProduct;
133+
dReverseProduct_dxi[b] *= reverseProduct;
134+
jacobian( a, b ) = -1/Keq * ( dReverseProduct_dxi[b] / forwardProduct - dForwardProduct_dxi[b] * reverseProduct / ( forwardProduct * forwardProduct ) );
130135
}
131136
else if constexpr( RESIDUAL_FORM == 2 )
132137
{
133-
jacobian( a, b ) = -dForwardProduct_dxi[b] / forwardProduct + dReverseProduct_dxi[b] / reverseProduct;
138+
// printf( "(%d) dForwardProduct_dxi[%d] = %f\n", a, b, dForwardProduct_dxi[b] );
139+
// printf( "(%d) dReverseProduct_dxi[%d] = %f\n", a, b, dReverseProduct_dxi[b] );
140+
jacobian( a, b ) = -dForwardProduct_dxi[b] + dReverseProduct_dxi[b];
141+
// printf( " jacobian( a, b ) = %g + %g = %g\n", -dForwardProduct_dxi[b], dReverseProduct_dxi[b], jacobian( a, b ) );
134142
}
135143
}
136144
}
@@ -162,7 +170,7 @@ EquilibriumReactions< REAL_TYPE,
162170
CArrayWrapper< double, numReactions, numReactions > jacobian;
163171

164172
REAL_TYPE residualNorm = 0.0;
165-
for( int k=0; k<100; ++k )
173+
for( int k=0; k<30; ++k )
166174
{
167175
computeResidualAndJacobian( temperature,
168176
params,
@@ -178,7 +186,7 @@ EquilibriumReactions< REAL_TYPE,
178186
}
179187
residualNorm = sqrt( residualNorm );
180188
printf( "iter, residualNorm = %2d, %16.10g \n", k, residualNorm );
181-
if( residualNorm < 1.0e-14 )
189+
if( residualNorm < 1.0e-12 )
182190
{
183191
printf( " converged\n" );
184192
break;
@@ -213,26 +221,28 @@ EquilibriumReactions< REAL_TYPE,
213221
printf( "},\n" );
214222
}
215223
printf( " }\n" );
216-
}
217224

218-
int numNonSymmetric = 0;
219-
for( int i=0; i<numReactions; ++i )
220-
{
221-
for( int j=i; j<numReactions; ++j )
225+
int numNonSymmetric = 0;
226+
for( int i=0; i<numReactions; ++i )
222227
{
223-
if( fabs( jacobian( i, j ) - jacobian( j, i ) ) > 0.5 * ( jacobian( i, j ) + jacobian( j, i ) )*1.0e-14 )
228+
for( int j=i; j<numReactions; ++j )
224229
{
225-
++numNonSymmetric;
226-
printf( "jacobian not symmetric: i = %d, j = %d, jacobian(i,j) = %g, jacobian(j,i) = %g\n", i, j, jacobian( i, j ), jacobian( j, i ) );
230+
if( fabs( jacobian( i, j ) - jacobian( j, i ) ) > 0.5 * ( jacobian( i, j ) + jacobian( j, i ) )*1.0e-14 )
231+
{
232+
++numNonSymmetric;
233+
printf( "jacobian not symmetric: i = %d, j = %d, jacobian(i,j) = %g, jacobian(j,i) = %g\n", i, j, jacobian( i, j ), jacobian( j, i ) );
234+
}
227235
}
228236
}
229-
}
230-
if( numNonSymmetric > 0 )
231-
{
232-
printf( "numNonSymmetric = %d\n", numNonSymmetric );
237+
if( numNonSymmetric > 0 )
238+
{
239+
printf( "numNonSymmetric = %d\n", numNonSymmetric );
240+
}
241+
242+
printf( " is J PD = %d\n", isPositiveDefinite< double, numReactions >( jacobian.data ) );
243+
233244
}
234245

235-
printf( " is J PD = %d\n", isPositiveDefinite< double, numReactions >( jacobian.data ) );
236246

237247
// solve for the change in xi
238248
solveNxN_Cholesky< double, numReactions >( jacobian.data, residual, dxi );

‎src/reactions/bulkGeneric/unitTests/testEquilibriumReactions.cpp‎

Lines changed: 19 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -150,12 +150,7 @@ void testEnforceEquilibrium( PARAMS_DATA const & params,
150150

151151
for( int r=0; r<numSpecies; ++r )
152152
{
153-
printf( "%18.12e\n", speciesConcentration[r] );
154-
}
155-
156-
for( int r=0; r<numSpecies; ++r )
157-
{
158-
EXPECT_NEAR( speciesConcentration[r], expectedSpeciesConcentrations[r], 1.0e-8 );
153+
EXPECT_NEAR( speciesConcentration[r], expectedSpeciesConcentrations[r], 1.0e-8 * expectedSpeciesConcentrations[r] );
159154
}
160155

161156
}
@@ -212,32 +207,26 @@ TEST( testEquilibriumReactions, testCarbonateSystem )
212207
};
213208

214209
double const expectedSpeciesConcentrations[18] =
215-
{ 9.884330823068e+06, // OH-
216-
4.217710930333e-04, // CO2
217-
2.114046643682e-01, // CO3-2
218-
1.640757112486e-01, // H2CO3
219-
1.341448126749e-08, // CaHCO3+
220-
1.953088434753e-07, // CaCO3
221-
3.224784300154e-07, // CaSO4
222-
1.447573069249e-02, // CaCl+
223-
2.253393830438e-02, // CaCl2
224-
1.692550572351e-06, // MgSO4
225-
4.520203771797e-03, // NaSO4-
226-
9.884331245975e+06, // H+
227-
9.764456680183e-05, // HCO3-
228-
1.689799801383e-03, // Ca+2
229-
2.757778119920e-02, // SO4-2
230-
1.830456392699e+00, // Cl-
231-
1.649830744943e-02, // Mg+2
232-
1.085479796228e+00 // Na+1
210+
{ 2.327841695586879e-11, // OH-
211+
3.745973700632716e-01, // CO2
212+
3.956656978189456e-11, // CO3-2
213+
9.629355924567627e-04, // H2CO3
214+
6.739226982791492e-05, // CaHCO3+
215+
1.065032288527957e-09, // CaCO3
216+
5.298329882666738e-03, // CaSO4
217+
5.844517547638333e-03, // CaCl+
218+
1.277319392670652e-02, // CaCl2
219+
6.618125707964991e-03, // MgSO4
220+
1.769217213462983e-02, // NaSO4-
221+
4.396954721488358e-04, // H+
222+
3.723009698453808e-04, // HCO3-
223+
1.471656530812871e-02, // Ca+2
224+
2.491372274738741e-03, // SO4-2
225+
1.858609094598949e+00, // Cl-
226+
9.881874292035110e-03, // Mg+2
227+
1.072307827865370e+00 // Na+1
233228
};
234229

235-
// std::cout<<" RESIDUAL_FORM 0:"<<std::endl;
236-
// testEnforceEquilibrium< double, 0 >( carbonateSystem.equilibriumReactionsParameters(),
237-
// initialSpeciesConcentration,
238-
// expectedSpeciesConcentrations );
239-
240-
241230
std::cout<<" RESIDUAL_FORM 2:"<<std::endl;
242231
testEnforceEquilibrium< double, 2 >( carbonateSystem.equilibriumReactionsParameters(),
243232
initialSpeciesConcentration,

0 commit comments

Comments
 (0)