Cantera  3.1.0
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 if (log_level == -7) {
219 log_level = 0;
220 } else {
221 warn_deprecated("MMCollisionInt::init", "The log_level parameter "
222 "is deprecated and will be removed after Cantera 3.1.");
223 }
224 m_loglevel = log_level;
225 debuglog("Collision Integral Polynomial Fits\n", m_loglevel > 0);
226 m_nmin = -1;
227 m_nmax = -1;
228
229 for (int n = 0; n < 37; n++) {
230 if (tsmin > tstar[n+1]) {
231 m_nmin = n;
232 }
233 if (tsmax > tstar[n+1]) {
234 m_nmax = n+1;
235 }
236 }
237 if (m_nmin < 0 || m_nmin >= 36 || m_nmax < 0 || m_nmax > 36) {
238 m_nmin = 0;
239 m_nmax = 36;
240 }
241 if (m_loglevel > 0) {
242 writelogf("T*_min = %g\n", tstar[m_nmin + 1]);
243 writelogf("T*_max = %g\n", tstar[m_nmax + 1]);
244 }
245 m_logTemp.resize(37);
246 double e22 = 0.0, ea = 0.0, eb = 0.0, ec = 0.0;
247
248 if (m_loglevel > 0) {
249 writelog("Collision integral fits at each tabulated T* vs. delta*.\n"
250 "These polynomial fits are used to interpolate between "
251 "columns (delta*)\n in the Monchick and Mason tables."
252 " They are only used for nonzero delta*.\n");
253 if (log_level < 4) {
254 writelog("polynomial coefficients not printed (log_level < 4)\n");
255 }
256 }
257
258 for (int i = 0; i < 37; i++) {
259 m_logTemp[i] = log(tstar[i+1]);
260 vector<double> c(DeltaDegree+1);
261
262 double rmserr = fitDelta(0, i, DeltaDegree, c.data());
263 if (log_level > 3) {
264 writelogf("\ndelta* fit at T* = %.6g\n", tstar[i+1]);
265 writelog("omega22 = [" + vec2str(c) + "]\n");
266 }
267 m_o22poly.push_back(c);
268 e22 = std::max(e22, rmserr);
269
270 rmserr = fitDelta(1, i, DeltaDegree, c.data());
271 m_apoly.push_back(c);
272 if (log_level > 3) {
273 writelog("A* = [" + vec2str(c) + "]\n");
274 }
275 ea = std::max(ea, rmserr);
276
277 rmserr = fitDelta(2, i, DeltaDegree, c.data());
278 m_bpoly.push_back(c);
279 if (log_level > 3) {
280 writelog("B* = [" + vec2str(c) + "]\n");
281 }
282 eb = std::max(eb, rmserr);
283
284 rmserr = fitDelta(3, i, DeltaDegree, c.data());
285 m_cpoly.push_back(c);
286 if (log_level > 3) {
287 writelog("C* = [" + vec2str(c) + "]\n");
288 }
289 ec = std::max(ec, rmserr);
290
291 if (log_level > 0) {
292 writelogf("max RMS errors in fits vs. delta*:\n"
293 " omega_22 = %12.6g \n"
294 " A* = %12.6g \n"
295 " B* = %12.6g \n"
296 " C* = %12.6g \n", e22, ea, eb, ec);
297 }
298 }
299}
300
301double MMCollisionInt::fitDelta(int table, int ntstar, int degree, double* c)
302{
303 vector<double> w(8);
304 double* begin = 0;
305 switch (table) {
306 case 0:
307 begin = omega22_table + 8*ntstar;
308 break;
309 case 1:
310 begin = astar_table + 8*(ntstar + 1);
311 break;
312 case 2:
313 begin = bstar_table + 8*(ntstar + 1);
314 break;
315 case 3:
316 begin = cstar_table + 8*(ntstar + 1);
317 break;
318 default:
319 return 0.0;
320 }
321 w[0] = -1.0;
322 return polyfit(8, degree, delta, begin, w.data(), c);
323}
324
325double MMCollisionInt::omega22(double ts, double deltastar)
326{
327 int i;
328 for (i = 0; i < 37; i++) {
329 if (ts < tstar22[i]) {
330 break;
331 }
332 }
333 int i1 = std::max(i - 1, 0);
334 int i2 = i1+3;
335 if (i2 > 36) {
336 i2 = 36;
337 i1 = i2 - 3;
338 }
339 vector<double> values(3);
340 for (i = i1; i < i2; i++) {
341 if (deltastar == 0.0) {
342 values[i-i1] = omega22_table[8*i];
343 } else {
344 values[i-i1] = poly6(deltastar, m_o22poly[i].data());
345 }
346 }
347 return quadInterp(log(ts), &m_logTemp[i1], values.data());
348}
349
350double MMCollisionInt::astar(double ts, double deltastar)
351{
352 int i;
353 for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
354 break;
355 }
356 int i1 = std::max(i - 1, 0);
357 int i2 = i1+3;
358 if (i2 > 36) {
359 i2 = 36;
360 i1 = i2 - 3;
361 }
362 vector<double> values(3);
363 for (i = i1; i < i2; i++) {
364 if (deltastar == 0.0) {
365 values[i-i1] = astar_table[8*(i + 1)];
366 } else {
367 values[i-i1] = poly6(deltastar, m_apoly[i].data());
368 }
369 }
370 return quadInterp(log(ts), &m_logTemp[i1], values.data());
371}
372
373double MMCollisionInt::bstar(double ts, double deltastar)
374{
375 int i;
376 for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
377 break;
378 }
379 int i1 = std::max(i - 1, 0);
380 int i2 = i1+3;
381 if (i2 > 36) {
382 i2 = 36;
383 i1 = i2 - 3;
384 }
385 vector<double> values(3);
386 for (i = i1; i < i2; i++) {
387 if (deltastar == 0.0) {
388 values[i-i1] = bstar_table[8*(i + 1)];
389 } else {
390 values[i-i1] = poly6(deltastar, m_bpoly[i].data());
391 }
392 }
393 return quadInterp(log(ts), &m_logTemp[i1], values.data());
394}
395
396double MMCollisionInt::cstar(double ts, double deltastar)
397{
398 int i;
399 for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
400 break;
401 }
402 int i1 = std::max(i - 1,0);
403 int i2 = i1+3;
404 if (i2 > 36) {
405 i2 = 36;
406 i1 = i2 - 3;
407 }
408 vector<double> values(3);
409 for (i = i1; i < i2; i++) {
410 if (deltastar == 0.0) {
411 values[i-i1] = cstar_table[8*(i + 1)];
412 } else {
413 values[i-i1] = poly6(deltastar, m_cpoly[i].data());
414 }
415 }
416 return quadInterp(log(ts), &m_logTemp[i1], values.data());
417}
418
419void MMCollisionInt::fit_omega22(int degree, double deltastar, double* o22)
420{
421 int n = m_nmax - m_nmin + 1;
422 vector<double> values(n);
423 vector<double> w(n);
424 double* logT = &m_logTemp[m_nmin];
425 for (int i = 0; i < n; i++) {
426 if (deltastar == 0.0) {
427 values[i] = omega22_table[8*(i + m_nmin)];
428 } else {
429 values[i] = poly6(deltastar, m_o22poly[i+m_nmin].data());
430 }
431 }
432 w[0]= -1.0;
433 double rmserr = polyfit(n, degree, logT, values.data(), w.data(), o22);
434 if (m_loglevel > 0 && rmserr > 0.01) {
435 warn_user("MMCollisionInt::fit_omega22",
436 "RMS error = {:12.6g} in omega_22 fit "
437 "with delta* = {:12.6g}", rmserr, deltastar);
438 }
439}
440
441void MMCollisionInt::fit(int degree, double deltastar, double* a, double* b, double* c)
442{
443 int n = m_nmax - m_nmin + 1;
444 vector<double> values(n);
445 vector<double> w(n);
446 double* logT = &m_logTemp[m_nmin];
447 for (int i = 0; i < n; i++) {
448 if (deltastar == 0.0) {
449 values[i] = astar_table[8*(i + m_nmin + 1)];
450 } else {
451 values[i] = poly6(deltastar, m_apoly[i+m_nmin].data());
452 }
453 }
454 w[0]= -1.0;
455 double rmserr = polyfit(n, degree, logT, values.data(), w.data(), a);
456
457 for (int i = 0; i < n; i++) {
458 if (deltastar == 0.0) {
459 values[i] = bstar_table[8*(i + m_nmin + 1)];
460 } else {
461 values[i] = poly6(deltastar, m_bpoly[i+m_nmin].data());
462 }
463 }
464 w[0]= -1.0;
465 rmserr = polyfit(n, degree, logT, values.data(), w.data(), b);
466
467 for (int i = 0; i < n; i++) {
468 if (deltastar == 0.0) {
469 values[i] = cstar_table[8*(i + m_nmin + 1)];
470 } else {
471 values[i] = poly6(deltastar, m_cpoly[i+m_nmin].data());
472 }
473 }
474 w[0]= -1.0;
475 rmserr = polyfit(n, degree, logT, values.data(), w.data(), c);
476 if (m_loglevel > 2) {
477 writelogf("\nT* fit at delta* = %.6g\n", deltastar);
478
479 writelog("astar = [" + vec2str(vector<double>(a, a+degree+1))+ "]\n");
480 if (rmserr > 0.01) {
481 warn_user("MMCollisionInt::fit",
482 "RMS error = {:12.6g} for A* fit", rmserr);
483 }
484
485 writelog("bstar = [" + vec2str(vector<double>(b, b+degree+1))+ "]\n");
486 if (rmserr > 0.01) {
487 warn_user("MMCollisionInt::fit",
488 "RMS error = {:12.6g} for B* fit", rmserr);
489 }
490
491 writelog("cstar = [" + vec2str(vector<double>(c, c+degree+1))+ "]\n");
492 if (rmserr > 0.01) {
493 warn_user("MMCollisionInt::fit",
494 "RMS error = {:12.6g} for C* fit", rmserr);
495 }
496 }
497}
498
499} // namespace
MMCollisionInt.h
Monchick and Mason collision integrals.
Cantera::MMCollisionInt::astar_table
static double astar_table[39 *8]
astar table from MM
Definition MMCollisionInt.h:69
Cantera::MMCollisionInt::cstar_table
static double cstar_table[39 *8]
cstar table from MM
Definition MMCollisionInt.h:75
Cantera::MMCollisionInt::tstar
static double tstar[39]
table of tstar values
Definition MMCollisionInt.h:66
Cantera::MMCollisionInt::omega22_table
static double omega22_table[37 *8]
Table of omega22 values from MM.
Definition MMCollisionInt.h:63
Cantera::MMCollisionInt::bstar_table
static double bstar_table[39 *8]
bstar table from MM
Definition MMCollisionInt.h:72
Cantera::MMCollisionInt::m_logTemp
vector< double > m_logTemp
Log temp.
Definition MMCollisionInt.h:78
Cantera::MMCollisionInt::init
void init(double tsmin, double tsmax, int loglevel=-7)
Initialize the object for calculation.
Definition MMCollisionInt.cpp:216
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
Cantera::warn_deprecated
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Definition AnyMap.cpp:1997
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...