HPCReact/src/reactions/bulkGeneric/EquilibriumReactionsAggregatePrimaryConcentration_impl.hpp at 0a2b5dc25f68fd213df51e041523cf03eb6d8054 · GEOS-DEV/HPCReact · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#if defined(__INTELLISENSE__)
#include "EquilibriumReactions.hpp"
#endif

#include "SpeciesUtilities.hpp"

namespace hpcReact
{
namespace bulkGeneric
{


template< typename REAL_TYPE,
          typename INT_TYPE,
          typename INDEX_TYPE >
template< typename PARAMS_DATA,
          typename ARRAY_1D,
          typename ARRAY_1D_TO_CONST,
          typename ARRAY_2D >
HPCREACT_HOST_DEVICE
inline
void
EquilibriumReactions< REAL_TYPE,
                      INT_TYPE,
                      INDEX_TYPE >::computeResidualAndJacobianAggregatePrimaryConcentrations( RealType const & temperature,
                                                                                              PARAMS_DATA const & params,
                                                                                              ARRAY_1D_TO_CONST const & targetAggregatePrimaryConcentrations,
                                                                                              ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration,
                                                                                              ARRAY_1D & residual,
                                                                                              ARRAY_2D & jacobian )
{
  HPCREACT_UNUSED_VAR( temperature );
  constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();

  RealType aggregatePrimaryConcentrations[numPrimarySpecies] = {0.0};
  ARRAY_2D dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations = {{{0.0}}};
  calculateAggregatePrimaryConcentrationsWrtLogC< REAL_TYPE, INT_TYPE, INDEX_TYPE >( params,
                                                                                     logPrimarySpeciesConcentration,
                                                                                     aggregatePrimaryConcentrations,
                                                                                     dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations );


  for( IndexType i=0; i<numPrimarySpecies; ++i )
  {
    residual[i] = -(1.0 - aggregatePrimaryConcentrations[i] / targetAggregatePrimaryConcentrations[i]);
    for( IndexType j=0; j<numPrimarySpecies; ++j )
    {
      jacobian( i, j ) = -dAggregatePrimarySpeciesConcentrationsDerivatives_dLogPrimarySpeciesConcentrations[i][j] / targetAggregatePrimaryConcentrations[i];
    }
  }
}

template< typename REAL_TYPE,
          typename INT_TYPE,
          typename INDEX_TYPE >
template< typename PARAMS_DATA,
          typename ARRAY_1D,
          typename ARRAY_1D_TO_CONST >
HPCREACT_HOST_DEVICE inline
void
EquilibriumReactions< REAL_TYPE,
                      INT_TYPE,
                      INDEX_TYPE >::enforceEquilibrium_Aggregate( REAL_TYPE const & temperature,
                                                                  PARAMS_DATA const & params,
                                                                  ARRAY_1D_TO_CONST const & logPrimarySpeciesConcentration0,
                                                                  ARRAY_1D & logPrimarySpeciesConcentration )
{
  HPCREACT_UNUSED_VAR( temperature );
  constexpr int numPrimarySpecies = PARAMS_DATA::numPrimarySpecies();

  double residual[numPrimarySpecies] = { 0.0 };
//  double aggregatePrimarySpeciesConcentration[numPrimarySpecies] = { 0.0 };
  double targetAggregatePrimarySpeciesConcentration[numPrimarySpecies] = { 0.0 };
  double dLogCp[numPrimarySpecies] = { 0.0 };
  CArrayWrapper< double, numPrimarySpecies, numPrimarySpecies > jacobian;


  for( int i=0; i<numPrimarySpecies; ++i )
  {
    targetAggregatePrimarySpeciesConcentration[i] = exp( logPrimarySpeciesConcentration0[i] );
    logPrimarySpeciesConcentration[i] = logPrimarySpeciesConcentration0[i];
  }


  REAL_TYPE residualNorm = 0.0;
  for( int k=0; k<30; ++k )
  {
    computeResidualAndJacobianAggregatePrimaryConcentrations( temperature,
                                                              params,
                                                              targetAggregatePrimarySpeciesConcentration,
                                                              logPrimarySpeciesConcentration,
                                                              residual,
                                                              jacobian );

    residualNorm = 0.0;
    for( int i = 0; i < numPrimarySpecies; ++i )
    {
      residualNorm += residual[i] * residual[i];
    }
    residualNorm = sqrt( residualNorm );
    printf( "iter, residualNorm = %2d, %16.10g \n", k, residualNorm );
    if( residualNorm < 1.0e-12 )
    {
      printf( " converged\n" );
      break;
    }

    solveNxN_pivoted< double, numPrimarySpecies >( jacobian.data, residual, dLogCp );


    for( IndexType i=0; i<numPrimarySpecies; ++i )
    {
      logPrimarySpeciesConcentration[i] += dLogCp[i];
    }

  }
}
} // namespace bulkGeneric
} // namespace hpcReact