Cantera  3.1.0b1
Loading...
Searching...
No Matches
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
14namespace Cantera
15{
16const int DeltaDegree = 6;
17
18double MMCollisionInt::delta[8] = {0.0, 0.25, 0.50, 0.75, 1.0,
19 1.5, 2.0, 2.5
20 };
21
22double 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
33double 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
40double 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
81double 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
88double 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
131double 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
174double 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
216void 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
295double 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
319double 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
344double 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
367double 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
390double 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
413void 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
435void 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:154
Cantera::writelogf
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:191
Cantera::writelog
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition global.h:171
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:263
Cantera
Namespace for the Cantera kernel.
Definition AnyMap.cpp:595
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...