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