Cantera  2.5.1
MMCollisionInt.cpp
Go to the documentation of this file.
1 /**
2  * @file MMCollisionInt.cpp
3  */
4 
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at https://cantera.org/license.txt for license and copyright information.
7 
8 #include "MMCollisionInt.h"
9 #include "cantera/base/utilities.h"
10 #include "cantera/numerics/polyfit.h"
11 #include "cantera/base/stringUtils.h"
12 
13 using namespace std;
14 
15 namespace Cantera
16 {
17 const int DeltaDegree = 6;
18 
19 double MMCollisionInt::delta[8] = {0.0, 0.25, 0.50, 0.75, 1.0,
20  1.5, 2.0, 2.5
21  };
22 
23 doublereal quadInterp(doublereal x0, doublereal* x, doublereal* y)
24 {
25  double dx21 = x[1] - x[0];
26  double dx32 = x[2] - x[1];
27  double dx31 = dx21 + dx32;
28  double dy32 = y[2] - y[1];
29  double dy21 = y[1] - y[0];
30  double a = (dx21*dy32 - dy21*dx32)/(dx21*dx31*dx32);
31  return a*(x0 - x[0])*(x0 - x[1]) + (dy21/dx21)*(x0 - x[1]) + y[1];
32 }
33 
34 double MMCollisionInt::tstar22[37] = {
35  0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
36  1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0,
37  5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0,
38  18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 50.0, 75.0, 100.0
39 };
40 
41 double MMCollisionInt::omega22_table[37*8] = {
42  4.1005, 4.266, 4.833, 5.742, 6.729, 8.624, 10.34, 11.89,
43  3.2626, 3.305, 3.516, 3.914, 4.433, 5.57, 6.637, 7.618,
44  2.8399, 2.836, 2.936, 3.168, 3.511, 4.329, 5.126, 5.874,
45  2.531, 2.522, 2.586, 2.749, 3.004, 3.64, 4.282, 4.895,
46  2.2837, 2.277, 2.329, 2.46, 2.665, 3.187, 3.727, 4.249,
47  2.0838, 2.081, 2.13, 2.243, 2.417, 2.862, 3.329, 3.786,
48  1.922, 1.924, 1.97, 2.072, 2.225, 2.614, 3.028, 3.435,
49  1.7902, 1.795, 1.84, 1.934, 2.07, 2.417, 2.788, 3.156,
50  1.6823, 1.689, 1.733, 1.82, 1.944, 2.258, 2.596, 2.933,
51  1.5929, 1.601, 1.644, 1.725, 1.838, 2.124, 2.435, 2.746,
52  1.4551, 1.465, 1.504, 1.574, 1.67, 1.913, 2.181, 2.451,
53  1.3551, 1.365, 1.4, 1.461, 1.544, 1.754, 1.989, 2.228,
54  1.28, 1.289, 1.321, 1.374, 1.447, 1.63, 1.838, 2.053,
55  1.2219, 1.231, 1.259, 1.306, 1.37, 1.532, 1.718, 1.912,
56  1.1757, 1.184, 1.209, 1.251, 1.307, 1.451, 1.618, 1.795,
57  1.0933, 1.1, 1.119, 1.15, 1.193, 1.304, 1.435, 1.578,
58  1.0388, 1.044, 1.059, 1.083, 1.117, 1.204, 1.31, 1.428,
59  0.99963, 1.004, 1.016, 1.035, 1.062, 1.133, 1.22, 1.319,
60  0.96988, 0.9732, 0.983, 0.9991, 1.021, 1.079, 1.153, 1.236,
61  0.92676, 0.9291, 0.936, 0.9473, 0.9628, 1.005, 1.058, 1.121,
62  0.89616, 0.8979, 0.903, 0.9114, 0.923, 0.9545, 0.9955, 1.044,
63  0.87272, 0.8741, 0.878, 0.8845, 0.8935, 0.9181, 0.9505, 0.9893,
64  0.85379, 0.8549, 0.858, 0.8632, 0.8703, 0.8901, 0.9164, 0.9482,
65  0.83795, 0.8388, 0.8414, 0.8456, 0.8515, 0.8678, 0.8895, 0.916,
66  0.82435, 0.8251, 0.8273, 0.8308, 0.8356, 0.8493, 0.8676, 0.8901,
67  0.80184, 0.8024, 0.8039, 0.8065, 0.8101, 0.8201, 0.8337, 0.8504,
68  0.78363, 0.784, 0.7852, 0.7872, 0.7899, 0.7976, 0.8081, 0.8212,
69  0.76834, 0.7687, 0.7696, 0.7712, 0.7733, 0.7794, 0.7878, 0.7983,
70  0.75518, 0.7554, 0.7562, 0.7575, 0.7592, 0.7642, 0.7711, 0.7797,
71  0.74364, 0.7438, 0.7445, 0.7455, 0.747, 0.7512, 0.7569, 0.7642,
72  0.71982, 0.72, 0.7204, 0.7211, 0.7221, 0.725, 0.7289, 0.7339,
73  0.70097, 0.7011, 0.7014, 0.7019, 0.7026, 0.7047, 0.7076, 0.7112,
74  0.68545, 0.6855, 0.6858, 0.6861, 0.6867, 0.6883, 0.6905, 0.6932,
75  0.67232, 0.6724, 0.6726, 0.6728, 0.6733, 0.6743, 0.6762, 0.6784,
76  0.65099, 0.651, 0.6512, 0.6513, 0.6516, 0.6524, 0.6534, 0.6546,
77  0.61397, 0.6141, 0.6143, 0.6145, 0.6147, 0.6148, 0.6148, 0.6147,
78  0.5887, 0.5889, 0.5894, 0.59, 0.5903, 0.5901, 0.5895, 0.5885
79 };
80 
81 // changed upper limit to 500 from 1.0e10 dgg 5/21/04
82 double MMCollisionInt::tstar[39] = {
83  0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
84  1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0,
85  5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0,
86  18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 50.0, 75.0, 100.0, 500.0
87 };
88 
89 double MMCollisionInt::astar_table[39*8] = {
90  1.0065, 1.0840, 1.0840, 1.0840, 1.0840, 1.0840, 1.0840, 1.0840,
91  1.0231, 1.0660, 1.0380, 1.0400, 1.0430, 1.0500, 1.0520, 1.0510,
92  1.0424, 1.0450, 1.0480, 1.0520, 1.0560, 1.0650, 1.0660, 1.0640,
93  1.0719, 1.0670, 1.0600, 1.0550, 1.0580, 1.0680, 1.0710, 1.0710,
94  1.0936, 1.0870, 1.0770, 1.0690, 1.0680, 1.0750, 1.0780, 1.0780,
95  1.1053, 1.0980, 1.0880, 1.0800, 1.0780, 1.0820, 1.0840, 1.0840,
96  1.1104, 1.1040, 1.0960, 1.0890, 1.0860, 1.0890, 1.0900, 1.0900,
97  1.1114, 1.1070, 1.1000, 1.0950, 1.0930, 1.0950, 1.0960, 1.0950,
98  1.1104, 1.1070, 1.1020, 1.0990, 1.0980, 1.1000, 1.1000, 1.0990,
99  1.1086, 1.1060, 1.1020, 1.1010, 1.1010, 1.1050, 1.1050, 1.1040,
100  1.1063, 1.1040, 1.1030, 1.1030, 1.1040, 1.1080, 1.1090, 1.1080,
101  1.1020, 1.1020, 1.1030, 1.1050, 1.1070, 1.1120, 1.1150, 1.1150,
102  1.0985, 1.0990, 1.1010, 1.1040, 1.1080, 1.1150, 1.1190, 1.1200,
103  1.0960, 1.0960, 1.0990, 1.1030, 1.1080, 1.1160, 1.1210, 1.1240,
104  1.0943, 1.0950, 1.0990, 1.1020, 1.1080, 1.1170, 1.1230, 1.1260,
105  1.0934, 1.0940, 1.0970, 1.1020, 1.1070, 1.1160, 1.1230, 1.1280,
106  1.0926, 1.0940, 1.0970, 1.0990, 1.1050, 1.1150, 1.1230, 1.1300,
107  1.0934, 1.0950, 1.0970, 1.0990, 1.1040, 1.1130, 1.1220, 1.1290,
108  1.0948, 1.0960, 1.0980, 1.1000, 1.1030, 1.1120, 1.1190, 1.1270,
109  1.0965, 1.0970, 1.0990, 1.1010, 1.1040, 1.1100, 1.1180, 1.1260,
110  1.0997, 1.1000, 1.1010, 1.1020, 1.1050, 1.1100, 1.1160, 1.1230,
111  1.1025, 1.1030, 1.1040, 1.1050, 1.1060, 1.1100, 1.1150, 1.1210,
112  1.1050, 1.1050, 1.1060, 1.1070, 1.1080, 1.1110, 1.1150, 1.1200,
113  1.1072, 1.1070, 1.1080, 1.1080, 1.1090, 1.1120, 1.1150, 1.1190,
114  1.1091, 1.1090, 1.1090, 1.1100, 1.1110, 1.1130, 1.1150, 1.1190,
115  1.1107, 1.1110, 1.1110, 1.1110, 1.1120, 1.1140, 1.1160, 1.1190,
116  1.1133, 1.1140, 1.1130, 1.1140, 1.1140, 1.1150, 1.1170, 1.1190,
117  1.1154, 1.1150, 1.1160, 1.1160, 1.1160, 1.1170, 1.1180, 1.1200,
118  1.1172, 1.1170, 1.1170, 1.1180, 1.1180, 1.1180, 1.1190, 1.1200,
119  1.1186, 1.1190, 1.1190, 1.1190, 1.1190, 1.1190, 1.1200, 1.1210,
120  1.1199, 1.1200, 1.1200, 1.1200, 1.1200, 1.1210, 1.1210, 1.1220,
121  1.1223, 1.1220, 1.1220, 1.1220, 1.1220, 1.1230, 1.1230, 1.1240,
122  1.1243, 1.1240, 1.1240, 1.1240, 1.1240, 1.1240, 1.1250, 1.1250,
123  1.1259, 1.1260, 1.1260, 1.1260, 1.1260, 1.1260, 1.1260, 1.1260,
124  1.1273, 1.1270, 1.1270, 1.1270, 1.1270, 1.1270, 1.1270, 1.1280,
125  1.1297, 1.1300, 1.1300, 1.1300, 1.1300, 1.1300, 1.1300, 1.1290,
126  1.1339, 1.1340, 1.1340, 1.1350, 1.1350, 1.1340, 1.1340, 1.1320,
127  1.1364, 1.1370, 1.1370, 1.1380, 1.1390, 1.1380, 1.1370, 1.1350,
128  1.14187, 1.14187, 1.14187, 1.14187, 1.14187, 1.14187, 1.14187,
129  1.14187
130 };
131 
132 double MMCollisionInt::bstar_table[39*8] = {
133  1.1852, 1.2963, 1.2963, 1.2963, 1.2963, 1.2963,1.2963, 1.2963,
134  1.1960, 1.216, 1.237, 1.269, 1.285, 1.290, 1.297, 1.294,
135  1.2451, 1.257, 1.340, 1.389, 1.366, 1.327, 1.314, 1.278,
136  1.2900, 1.294, 1.272, 1.258, 1.262, 1.282, 1.290, 1.299,
137  1.2986, 1.291, 1.284, 1.278, 1.277, 1.288, 1.294, 1.297,
138  1.2865, 1.281, 1.276, 1.272, 1.277, 1.286, 1.292, 1.298,
139  1.2665, 1.264, 1.261, 1.263, 1.269, 1.284, 1.292, 1.298,
140  1.2455, 1.244, 1.248, 1.255, 1.262, 1.278, 1.289, 1.296,
141  1.2253, 1.225, 1.234, 1.240, 1.252, 1.271, 1.284, 1.295,
142  1.2078, 1.210, 1.216, 1.227, 1.242, 1.264, 1.281, 1.292,
143  1.1919, 1.192, 1.205, 1.216, 1.230, 1.256, 1.273, 1.287,
144  1.1678, 1.172, 1.181, 1.195, 1.209, 1.237, 1.261, 1.277,
145  1.1496, 1.155, 1.161, 1.174, 1.189, 1.221, 1.246, 1.266,
146  1.1366, 1.141, 1.147, 1.159, 1.174, 1.202, 1.231, 1.256,
147  1.1270, 1.130, 1.138, 1.148, 1.162, 1.191, 1.218, 1.242,
148  1.1197, 1.122, 1.129, 1.140, 1.149, 1.178, 1.205, 1.231,
149  1.1080, 1.110, 1.116, 1.122, 1.132, 1.154, 1.180, 1.205,
150  1.1016, 1.103, 1.107, 1.112, 1.120, 1.138, 1.160, 1.183,
151  1.0980, 1.099, 1.102, 1.106, 1.112, 1.127, 1.145, 1.165,
152  1.0958, 1.097, 1.099, 1.102, 1.107, 1.119, 1.135, 1.153,
153  1.0935, 1.094, 1.095, 1.097, 1.100, 1.109, 1.120, 1.134,
154  1.0925, 1.092, 1.094, 1.095, 1.098, 1.104, 1.112, 1.122,
155  1.0922, 1.092, 1.093, 1.094, 1.096, 1.100, 1.106, 1.115,
156  1.0922, 1.092, 1.093, 1.093, 1.095, 1.098, 1.103, 1.110,
157  1.0923, 1.092, 1.093, 1.093, 1.094, 1.097, 1.101, 1.106,
158  1.0923, 1.092, 1.092, 1.093, 1.094, 1.096, 1.099, 1.103,
159  1.0927, 1.093, 1.093, 1.093, 1.094, 1.095, 1.098, 1.101,
160  1.0930, 1.093, 1.093, 1.093, 1.094, 1.094, 1.096, 1.099,
161  1.0933, 1.094, 1.093, 1.094, 1.094, 1.095, 1.096, 1.098,
162  1.0937, 1.093, 1.094, 1.094, 1.094, 1.094, 1.096, 1.097,
163  1.0939, 1.094, 1.094, 1.094, 1.094, 1.095, 1.095, 1.097,
164  1.0943, 1.094, 1.094, 1.094, 1.095, 1.095, 1.096, 1.096,
165  1.0944, 1.095, 1.094, 1.094, 1.094, 1.095, 1.095, 1.096,
166  1.0944, 1.094, 1.095, 1.094, 1.094, 1.095, 1.096, 1.096,
167  1.0943, 1.095, 1.094, 1.094, 1.095, 1.095, 1.095, 1.095,
168  1.0941, 1.094, 1.094, 1.094, 1.094, 1.094, 1.094, 1.096,
169  1.0947, 1.095, 1.094, 1.094, 1.093, 1.093, 1.094, 1.095,
170  1.0957, 1.095, 1.094, 1.093, 1.092, 1.093, 1.093, 1.094,
171  1.10185, 1.10185, 1.10185, 1.10185, 1.10185, 1.10185, 1.10185,
172  1.10185
173 };
174 
175 double MMCollisionInt::cstar_table[39*8] = {
176  0.8889, 0.77778, 0.77778,0.77778,0.77778,0.77778,0.77778,0.77778,
177  0.88575, 0.8988, 0.8378, 0.8029, 0.7876, 0.7805, 0.7799, 0.7801,
178  0.87268, 0.8692,0.8647,0.8479,0.8237,0.7975,0.7881,0.7784,
179  0.85182, 0.8525,0.8366,0.8198,0.8054,0.7903,0.7839,0.782,
180  0.83542, 0.8362,0.8306,0.8196,0.8076,0.7918,0.7842,0.7806,
181  0.82629, 0.8278,0.8252,0.8169,0.8074,0.7916,0.7838,0.7802,
182  0.82299, 0.8249,0.823,0.8165,0.8072,0.7922,0.7839,0.7798,
183  0.82357, 0.8257,0.8241,0.8178,0.8084,0.7927,0.7839,0.7794,
184  0.82657, 0.828,0.8264,0.8199,0.8107,0.7939,0.7842,0.7796,
185  0.8311, 0.8234,0.8295,0.8228,0.8136,0.796,0.7854,0.7798,
186  0.8363, 0.8366,0.8342,0.8267,0.8168,0.7986,0.7864,0.7805,
187  0.84762, 0.8474,0.8438,0.8358,0.825,0.8041,0.7904,0.7822,
188  0.85846, 0.8583,0.853,0.8444,0.8336,0.8118,0.7957,0.7854,
189  0.8684, 0.8674,0.8619,0.8531,0.8423,0.8186,0.8011,0.7898,
190  0.87713, 0.8755,0.8709,0.8616,0.8504,0.8265,0.8072,0.7939,
191  0.88479, 0.8831,0.8779,0.8695,0.8578,0.8338,0.8133,0.799,
192  0.89972, 0.8986,0.8936,0.8846,0.8742,0.8504,0.8294,0.8125,
193  0.91028, 0.9089,0.9043,0.8967,0.8869,0.8649,0.8438,0.8253,
194  0.91793, 0.9166,0.9125,0.9058,0.897,0.8768,0.8557,0.8372,
195  0.92371, 0.9226,0.9189,0.9128,0.905,0.8861,0.8664,0.8484,
196  0.93135, 0.9304,0.9274,0.9226,0.9164,0.9006,0.8833,0.8662,
197  0.93607, 0.9353,0.9329,0.9291,0.924,0.9109,0.8958,0.8802,
198  0.93927, 0.9387,0.9366,0.9334,0.9292,0.9162,0.905,0.8911,
199  0.94149, 0.9409,0.9393,0.9366,0.9331,0.9236,0.9122,0.8997,
200  0.94306, 0.9426,0.9412,0.9388,0.9357,0.9276,0.9175,0.9065,
201  0.94419, 0.9437,0.9425,0.9406,0.938,0.9308,0.9219,0.9119,
202  0.94571, 0.9455,0.9445,0.943,0.9409,0.9353,0.9283,0.9201,
203  0.94662, 0.9464,0.9456,0.9444,0.9428,0.9382,0.9325,0.9258,
204  0.94723, 0.9471,0.9464,0.9455,0.9442,0.9405,0.9355,0.9298,
205  0.94764, 0.9474,0.9469,0.9462,0.945,0.9418,0.9378,0.9328,
206  0.9479, 0.9478,0.9474,0.9465,0.9457,0.943,0.9394,0.9352,
207  0.94827, 0.9481,0.948,0.9472,0.9467,0.9447,0.9422,0.9391,
208  0.94842, 0.9484,0.9481,0.9478,0.9472,0.9458,0.9437,0.9415,
209  0.94852, 0.9484,0.9483,0.948,0.9475,0.9465,0.9449,0.943,
210  0.94861, 0.9487,0.9484,0.9481,0.9479,0.9468,0.9455,0.943,
211  0.94872, 0.9486,0.9486,0.9483,0.9482,0.9475,0.9464,0.9452,
212  0.94881, 0.9488,0.9489,0.949,0.9487,0.9482,0.9476,0.9468,
213  0.94863, 0.9487,0.9489,0.9491,0.9493,0.9491,0.9483,0.9476,
214  0.94444, 0.94444,0.94444,0.94444,0.94444,0.94444,0.94444,0.94444
215 };
216 
217 void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level)
218 {
219  m_loglevel = log_level;
220  debuglog("Collision Integral Polynomial Fits\n", m_loglevel > 0);
221  m_nmin = -1;
222  m_nmax = -1;
223 
224  for (int n = 0; n < 37; n++) {
225  if (tsmin > tstar[n+1]) {
226  m_nmin = n;
227  }
228  if (tsmax > tstar[n+1]) {
229  m_nmax = n+1;
230  }
231  }
232  if (m_nmin < 0 || m_nmin >= 36 || m_nmax < 0 || m_nmax > 36) {
233  m_nmin = 0;
234  m_nmax = 36;
235  }
236  if (m_loglevel > 0) {
237  writelogf("T*_min = %g\n", tstar[m_nmin + 1]);
238  writelogf("T*_max = %g\n", tstar[m_nmax + 1]);
239  }
240  m_logTemp.resize(37);
241  doublereal e22 = 0.0, ea = 0.0, eb = 0.0, ec = 0.0;
242 
243  if (m_loglevel > 0) {
244  writelog("Collision integral fits at each tabulated T* vs. delta*.\n"
245  "These polynomial fits are used to interpolate between "
246  "columns (delta*)\n in the Monchick and Mason tables."
247  " They are only used for nonzero delta*.\n");
248  if (log_level < 4) {
249  writelog("polynomial coefficients not printed (log_level < 4)\n");
250  }
251  }
252 
253  for (int i = 0; i < 37; i++) {
254  m_logTemp[i] = log(tstar[i+1]);
255  vector_fp c(DeltaDegree+1);
256 
257  double rmserr = fitDelta(0, i, DeltaDegree, c.data());
258  if (log_level > 3) {
259  writelogf("\ndelta* fit at T* = %.6g\n", tstar[i+1]);
260  writelog("omega22 = [" + vec2str(c) + "]\n");
261  }
262  m_o22poly.push_back(c);
263  e22 = std::max(e22, rmserr);
264 
265  rmserr = fitDelta(1, i, DeltaDegree, c.data());
266  m_apoly.push_back(c);
267  if (log_level > 3) {
268  writelog("A* = [" + vec2str(c) + "]\n");
269  }
270  ea = std::max(ea, rmserr);
271 
272  rmserr = fitDelta(2, i, DeltaDegree, c.data());
273  m_bpoly.push_back(c);
274  if (log_level > 3) {
275  writelog("B* = [" + vec2str(c) + "]\n");
276  }
277  eb = std::max(eb, rmserr);
278 
279  rmserr = fitDelta(3, i, DeltaDegree, c.data());
280  m_cpoly.push_back(c);
281  if (log_level > 3) {
282  writelog("C* = [" + vec2str(c) + "]\n");
283  }
284  ec = std::max(ec, rmserr);
285 
286  if (log_level > 0) {
287  writelogf("max RMS errors in fits vs. delta*:\n"
288  " omega_22 = %12.6g \n"
289  " A* = %12.6g \n"
290  " B* = %12.6g \n"
291  " C* = %12.6g \n", e22, ea, eb, ec);
292  }
293  }
294 }
295 
296 doublereal MMCollisionInt::fitDelta(int table, int ntstar, int degree, doublereal* c)
297 {
298  vector_fp w(8);
299  doublereal* begin = 0;
300  switch (table) {
301  case 0:
302  begin = omega22_table + 8*ntstar;
303  break;
304  case 1:
305  begin = astar_table + 8*(ntstar + 1);
306  break;
307  case 2:
308  begin = bstar_table + 8*(ntstar + 1);
309  break;
310  case 3:
311  begin = cstar_table + 8*(ntstar + 1);
312  break;
313  default:
314  return 0.0;
315  }
316  w[0] = -1.0;
317  return polyfit(8, degree, delta, begin, w.data(), c);
318 }
319 
320 doublereal MMCollisionInt::omega22(double ts, double deltastar)
321 {
322  int i;
323  for (i = 0; i < 37; i++) {
324  if (ts < tstar22[i]) {
325  break;
326  }
327  }
328  int i1 = std::max(i - 1, 0);
329  int i2 = i1+3;
330  if (i2 > 36) {
331  i2 = 36;
332  i1 = i2 - 3;
333  }
334  vector_fp values(3);
335  for (i = i1; i < i2; i++) {
336  if (deltastar == 0.0) {
337  values[i-i1] = omega22_table[8*i];
338  } else {
339  values[i-i1] = poly5(deltastar, m_o22poly[i].data());
340  }
341  }
342  return quadInterp(log(ts), &m_logTemp[i1], values.data());
343 }
344 
345 doublereal MMCollisionInt::astar(double ts, double deltastar)
346 {
347  int i;
348  for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
349  break;
350  }
351  int i1 = std::max(i - 1, 0);
352  int i2 = i1+3;
353  if (i2 > 36) {
354  i2 = 36;
355  i1 = i2 - 3;
356  }
357  vector_fp values(3);
358  for (i = i1; i < i2; i++) {
359  if (deltastar == 0.0) {
360  values[i-i1] = astar_table[8*(i + 1)];
361  } else {
362  values[i-i1] = poly5(deltastar, m_apoly[i].data());
363  }
364  }
365  return quadInterp(log(ts), &m_logTemp[i1], values.data());
366 }
367 
368 doublereal MMCollisionInt::bstar(double ts, double deltastar)
369 {
370  int i;
371  for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
372  break;
373  }
374  int i1 = std::max(i - 1, 0);
375  int i2 = i1+3;
376  if (i2 > 36) {
377  i2 = 36;
378  i1 = i2 - 3;
379  }
380  vector_fp values(3);
381  for (i = i1; i < i2; i++) {
382  if (deltastar == 0.0) {
383  values[i-i1] = bstar_table[8*(i + 1)];
384  } else {
385  values[i-i1] = poly5(deltastar, m_bpoly[i].data());
386  }
387  }
388  return quadInterp(log(ts), &m_logTemp[i1], values.data());
389 }
390 
391 doublereal MMCollisionInt::cstar(double ts, double deltastar)
392 {
393  int i;
394  for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
395  break;
396  }
397  int i1 = std::max(i - 1,0);
398  int i2 = i1+3;
399  if (i2 > 36) {
400  i2 = 36;
401  i1 = i2 - 3;
402  }
403  vector_fp values(3);
404  for (i = i1; i < i2; i++) {
405  if (deltastar == 0.0) {
406  values[i-i1] = cstar_table[8*(i + 1)];
407  } else {
408  values[i-i1] = poly5(deltastar, m_cpoly[i].data());
409  }
410  }
411  return quadInterp(log(ts), &m_logTemp[i1], values.data());
412 }
413 
414 void MMCollisionInt::fit_omega22(int degree, doublereal deltastar,
415  doublereal* o22)
416 {
417  int n = m_nmax - m_nmin + 1;
418  vector_fp values(n);
419  vector_fp w(n);
420  doublereal* logT = &m_logTemp[m_nmin];
421  for (int i = 0; i < n; i++) {
422  if (deltastar == 0.0) {
423  values[i] = omega22_table[8*(i + m_nmin)];
424  } else {
425  values[i] = poly5(deltastar, m_o22poly[i+m_nmin].data());
426  }
427  }
428  w[0]= -1.0;
429  double rmserr = polyfit(n, degree, logT, values.data(), w.data(), o22);
430  if (m_loglevel > 0 && rmserr > 0.01) {
431  warn_user("MMCollisionInt::fit_omega22",
432  "RMS error = {:12.6g} in omega_22 fit "
433  "with delta* = {:12.6g}", rmserr, deltastar);
434  }
435 }
436 
437 void MMCollisionInt::fit(int degree, doublereal deltastar,
438  doublereal* a, doublereal* b, doublereal* c)
439 {
440  int n = m_nmax - m_nmin + 1;
441  vector_fp values(n);
442  vector_fp w(n);
443  doublereal* logT = &m_logTemp[m_nmin];
444  for (int i = 0; i < n; i++) {
445  if (deltastar == 0.0) {
446  values[i] = astar_table[8*(i + m_nmin + 1)];
447  } else {
448  values[i] = poly5(deltastar, m_apoly[i+m_nmin].data());
449  }
450  }
451  w[0]= -1.0;
452  double rmserr = polyfit(n, degree, logT, values.data(), w.data(), a);
453 
454  for (int i = 0; i < n; i++) {
455  if (deltastar == 0.0) {
456  values[i] = bstar_table[8*(i + m_nmin + 1)];
457  } else {
458  values[i] = poly5(deltastar, m_bpoly[i+m_nmin].data());
459  }
460  }
461  w[0]= -1.0;
462  rmserr = polyfit(n, degree, logT, values.data(), w.data(), b);
463 
464  for (int i = 0; i < n; i++) {
465  if (deltastar == 0.0) {
466  values[i] = cstar_table[8*(i + m_nmin + 1)];
467  } else {
468  values[i] = poly5(deltastar, m_cpoly[i+m_nmin].data());
469  }
470  }
471  w[0]= -1.0;
472  rmserr = polyfit(n, degree, logT, values.data(), w.data(), c);
473  if (m_loglevel > 2) {
474  writelogf("\nT* fit at delta* = %.6g\n", deltastar);
475 
476  writelog("astar = [" + vec2str(vector_fp(a, a+degree+1))+ "]\n");
477  if (rmserr > 0.01) {
478  warn_user("MMCollisionInt::fit",
479  "RMS error = {:12.6g} for A* fit", rmserr);
480  }
481 
482  writelog("bstar = [" + vec2str(vector_fp(b, b+degree+1))+ "]\n");
483  if (rmserr > 0.01) {
484  warn_user("MMCollisionInt::fit",
485  "RMS error = {:12.6g} for B* fit", rmserr);
486  }
487 
488  writelog("cstar = [" + vec2str(vector_fp(c, c+degree+1))+ "]\n");
489  if (rmserr > 0.01) {
490  warn_user("MMCollisionInt::fit",
491  "RMS error = {:12.6g} for C* fit", rmserr);
492  }
493  }
494 }
495 
496 } // namespace
MMCollisionInt.h
Monk and Monchick collision integrals.
Cantera::debuglog
void debuglog(const std::string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition: global.h:140
Cantera::warn_user
void warn_user(const std::string &method, const std::string &msg, const Args &... args)
Definition: global.h:206
Cantera::vector_fp
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:180
Cantera::writelog
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:158
Cantera::writelogf
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:179
Cantera
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
Cantera::vec2str
std::string vec2str(const vector_fp &v, const std::string &fmt, const std::string &sep)
Convert a vector to a string (separated by commas)
Definition: stringUtils.cpp:34
Cantera::poly5
R poly5(D x, R *c)
Templated evaluation of a polynomial of order 5.
Definition: utilities.h:479
Cantera::polyfit
double polyfit(size_t n, size_t deg, const double *xp, const double *yp, const double *wp, double *pp)
Fits a polynomial function to a set of data points.
Definition: polyfit.cpp:14
stringUtils.h
Contains declarations for string manipulation functions within Cantera.
utilities.h
Various templated functions that carry out common vector operations (see Templated Utility Functions)...
Generated by doxygen 1.9.1