1 #if !defined(NO_TABULAR_BACKENDS)
9 const std::vector<std::vector<CellCoeffs>>& coeffs,
double x,
double y,
10 std::size_t& i, std::size_t& j) {
19 throw ValueError(
format(
"Cell is invalid and has no good neighbors for x = %g, y= %g", x, y));
28 const double otherval, std::size_t& i, std::size_t& j) {
37 throw ValueError(
format(
"Cell is invalid and has no good neighbors for x = %g, y = %g", value1, otherval));
51 std::size_t i, std::size_t j) {
53 std::vector<std::vector<double>>* f = NULL;
62 throw ValueError(
format(
"invalid output variable to BicubicBackend::evaluate_single_phase_transport"));
64 double x1 = table.
xvec[i], x2 = table.
xvec[i + 1], y1 = table.
yvec[j], y2 = table.
yvec[j + 1];
65 double f11 = (*f)[i][j], f12 = (*f)[i][j + 1], f21 = (*f)[i + 1][j], f22 = (*f)[i + 1][j + 1];
67 1 / ((x2 - x1) * (y2 - y1)) * (f11 * (x2 - x) * (y2 - y) + f21 * (x - x1) * (y2 - y) + f12 * (x2 - x) * (y - y1) + f22 * (x - x1) * (y - y1));
78 throw ValueError(
"Invalid output variable in evaluate_single_phase_transport");
84 const parameters output,
const double x,
const double y,
const std::size_t i,
85 const std::size_t j) {
90 const std::vector<double>& alpha = cell.
get(output);
93 double xhat = (x - table.
xvec[i]) / (table.
xvec[i + 1] - table.
xvec[i]);
94 double yhat = (y - table.
yvec[j]) / (table.
yvec[j + 1] - table.
yvec[j]);
98 double B0 = ((((0) + alpha[3 * 4 + 0]) * yhat + alpha[2 * 4 + 0]) * yhat + alpha[1 * 4 + 0]) * yhat + alpha[0 * 4 + 0];
100 double B1 = ((((0) + alpha[3 * 4 + 1]) * yhat + alpha[2 * 4 + 1]) * yhat + alpha[1 * 4 + 1]) * yhat + alpha[0 * 4 + 1];
102 double B2 = ((((0) + alpha[3 * 4 + 2]) * yhat + alpha[2 * 4 + 2]) * yhat + alpha[1 * 4 + 2]) * yhat + alpha[0 * 4 + 2];
104 double B3 = ((((0) + alpha[3 * 4 + 3]) * yhat + alpha[2 * 4 + 3]) * yhat + alpha[1 * 4 + 3]) * yhat + alpha[0 * 4 + 3];
106 double val = ((((0) + B3) * xhat + B2) * xhat + B1) * xhat + B0;
126 throw ValueError(
"Invalid output variable in evaluate_single_phase");
132 parameters output,
double x,
double y, std::size_t i, std::size_t j, std::size_t Nx,
139 const std::vector<double>& alpha = cell.
get(output);
142 double xhat = (x - table.
xvec[i]) / (table.
xvec[i + 1] - table.
xvec[i]);
143 double yhat = (y - table.
yvec[j]) / (table.
yvec[j + 1] - table.
yvec[j]);
144 double dxhatdx = 1 / (table.
xvec[i + 1] - table.
xvec[i]);
145 double dyhatdy = 1 / (table.
yvec[j + 1] - table.
yvec[j]);
149 if (Nx == 1 && Ny == 0) {
150 if (output == table.
xkey) {
153 if (output == table.
ykey) {
156 for (std::size_t l = 1; l < 4; ++l) {
157 for (std::size_t m = 0; m < 4; ++m) {
158 val += alpha[m * 4 + l] * l * pow(xhat,
static_cast<int>(l - 1)) * pow(yhat,
static_cast<int>(m));
162 return val * dxhatdx;
163 }
else if (Ny == 1 && Nx == 0) {
164 if (output == table.
ykey) {
167 if (output == table.
xkey) {
170 for (std::size_t l = 0; l < 4; ++l) {
171 for (std::size_t m = 1; m < 4; ++m) {
172 val += alpha[m * 4 + l] * pow(xhat,
static_cast<int>(l)) * m * pow(yhat,
static_cast<int>(m - 1));
176 return val * dyhatdy;
184 parameters other_key,
double other,
double y, std::size_t i, std::size_t j) {
189 const std::vector<double>& alpha = cell.
get(other_key);
192 double yhat = (y - table.
yvec[j]) / (table.
yvec[j + 1] - table.
yvec[j]);
194 double y_0 = 1, y_1 = yhat, y_2 = yhat * yhat, y_3 = yhat * yhat * yhat;
196 double a = alpha[3 + 0 * 4] * y_0 + alpha[3 + 1 * 4] * y_1 + alpha[3 + 2 * 4] * y_2 + alpha[3 + 3 * 4] * y_3;
197 double b = alpha[2 + 0 * 4] * y_0 + alpha[2 + 1 * 4] * y_1 + alpha[2 + 2 * 4] * y_2 + alpha[2 + 3 * 4] * y_3;
198 double c = alpha[1 + 0 * 4] * y_0 + alpha[1 + 1 * 4] * y_1 + alpha[1 + 2 * 4] * y_2 + alpha[1 + 3 * 4] * y_3;
199 double d = alpha[0 + 0 * 4] * y_0 + alpha[0 + 1 * 4] * y_1 + alpha[0 + 2 * 4] * y_2 + alpha[0 + 3 * 4] * y_3 - other;
201 double xhat0, xhat1, xhat2, val, xhat = _HUGE;
206 xhat = std::abs(xhat0) < std::abs(xhat1) ? xhat0 : xhat1;
208 if (std::abs(xhat0) < std::abs(xhat1) && std::abs(xhat0) < std::abs(xhat2)) {
212 else if (std::abs(xhat1) < std::abs(xhat2)) {
218 throw ValueError(
"Could not find a solution in invert_single_phase_x");
223 val = xhat * (table.
xvec[i + 1] - table.
xvec[i]) + table.
xvec[i];
226 switch (table.
xkey) {
234 throw ValueError(
"Invalid output variable in invert_single_phase_x");
240 parameters other_key,
double other,
double x, std::size_t i, std::size_t j) {
245 const std::vector<double>& alpha = cell.
get(other_key);
248 double xhat = (x - table.
xvec[i]) / (table.
xvec[i + 1] - table.
xvec[i]);
250 double x_0 = 1, x_1 = xhat, x_2 = xhat * xhat, x_3 = xhat * xhat * xhat;
252 double a = alpha[0 + 3 * 4] * x_0 + alpha[1 + 3 * 4] * x_1 + alpha[2 + 3 * 4] * x_2 + alpha[3 + 3 * 4] * x_3;
253 double b = alpha[0 + 2 * 4] * x_0 + alpha[1 + 2 * 4] * x_1 + alpha[2 + 2 * 4] * x_2 + alpha[3 + 2 * 4] * x_3;
254 double c = alpha[0 + 1 * 4] * x_0 + alpha[1 + 1 * 4] * x_1 + alpha[2 + 1 * 4] * x_2 + alpha[3 + 1 * 4] * x_3;
255 double d = alpha[0 + 0 * 4] * x_0 + alpha[1 + 0 * 4] * x_1 + alpha[2 + 0 * 4] * x_2 + alpha[3 + 0 * 4] * x_3 - other;
257 double yhat0, yhat1, yhat2, val, yhat = _HUGE;
262 yhat = std::abs(yhat0) < std::abs(yhat1) ? yhat0 : yhat1;
264 if (std::abs(yhat0) < std::abs(yhat1) && std::abs(yhat0) < std::abs(yhat2)) {
268 else if (std::abs(yhat1) < std::abs(yhat2)) {
274 throw ValueError(
"Could not find a solution in invert_single_phase_x");
279 val = yhat * (table.
yvec[j + 1] - table.
yvec[j]) + table.
yvec[j];
282 switch (table.
ykey) {
287 throw ValueError(
"Invalid output variable in invert_single_phase_x");