From f03aeabc388fbdfea46ff3d95fbb09dcb6049dd3 Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Fri, 20 Dec 2019 14:24:31 +0100 Subject: [PATCH] Refactored cell volume calculations for speed --- .../utility/numeric/calculateCellVol.cpp | 156 ++++++++---------- 1 file changed, 69 insertions(+), 87 deletions(-) diff --git a/src/opm/common/utility/numeric/calculateCellVol.cpp b/src/opm/common/utility/numeric/calculateCellVol.cpp index daece892a..51da54137 100755 --- a/src/opm/common/utility/numeric/calculateCellVol.cpp +++ b/src/opm/common/utility/numeric/calculateCellVol.cpp @@ -40,112 +40,94 @@ */ double C(const double* r, int i1, int i2, int i3){ + int g = i1 + i2 * 2 + i3 * 4; - double temp; + if (g == 0) + return r[0]; - temp = 0.0; + if (g == 1) + return r[1] - r[0]; - if ((i1 == 0) && (i2 == 0) && (i3 == 0)){ - temp = r[0]; - } + if (g == 2) + return r[2] - r[0]; - if ((i1 == 1) && (i2 == 0) && (i3 == 0)){ - temp = r[1] - r[0]; - } + if (g == 3) + return r[3] + r[0] - r[2] - r[1]; - if ((i1 == 0) && (i2 == 1) && (i3 == 0)){ - temp = r[2] - r[0]; - } + if (g == 4) + return r[4] - r[0]; - if ((i1 == 0) && (i2 == 0) && (i3 == 1)){ - temp = r[4]-r[0]; - } + if (g == 5) + return r[5] + r[0] - r[4] - r[1]; - if ((i1 == 1) && (i2 == 1) && (i3 == 0)){ - temp = r[3] + r[0] - r[2] - r[1]; - } + if (g == 6) + return r[6] + r[0] - r[4] - r[2]; - if ((i1 == 0) && (i2 == 1) && (i3 == 1)){ - temp = r[6] + r[0] - r[4] - r[2]; - } - - if ((i1 == 1) && (i2 == 0) && (i3 == 1)){ - temp = r[5] + r[0] - r[4] - r[1]; - } - - if ((i1 == 1) && (i2 == 1) && (i3 == 1)){ - temp = r[7] + r[4] + r[2] + r[1] - r[6] - r[5] - r[3] - r[0]; - } - - return temp; + return r[7] + r[4] + r[2] + r[1] - r[6] - r[5] - r[3] - r[0]; } -double perm123sign(int i1, int i2, int i3){ +struct pqr_t { + int pb; + int pg; + int qa; + int qg; + int ra; + int rb; +}; - double temp; - - if ((i1 ==1 ) && (i2==2) && (i3 == 3)){ - temp = 1.0; - } else if ((i1 == 1) && (i2==3) && (i3 == 2)){ - temp = -1.0; - } else if ((i1 == 2) && (i2 == 1) && (i3 == 3)){ - temp = -1.0; - } else if ((i1 == 2) && (i2 == 3) && (i3 == 1)){ - temp = 1.0; - } else if ((i1 == 3) && (i2 == 1) && (i3 == 2)){ - temp = 1.0; - } else if ((i1 == 3) && (i2 == 2) && (i3 == 1)){ - temp = -1.0; - } else { - OPM_THROW(std::logic_error, "Wrong indices in perm123sign"); - } - - return temp; -} double calculateCellVol(const std::array& X, const std::array& Y, const std::array& Z){ + /* + The permutation array should be ordered so that the sign: - double volume = 0.0; - const double* vect[3]; - int permutation[] = {1, 2, 3}; + sign = (-1)^N, N = # permutations - do { - for (int j = 0; j < 3; ++j){ - if (permutation[j] == 1){ - vect[j] = X.data(); - } else if (permutation[j] == 2){ - vect[j] = Y.data(); - } else if (permutation[j] == 3){ - vect[j] = Z.data(); - } else { - // this condition can never happen, since all values in 'permutation' - // is covered, but compiler analysis may not be deep enough to not give - // warnings about 'vect' being uninitialized further down. - assert(false); - vect[j] = 0; - } - } + is alternating - so that the sign can just be changed multiplying with -1. + */ + static const std::array< std::array, 6 > permutation = {{{ 0, 1, 2}, + { 0, 2, 1}, + { 1, 2, 0}, + { 1, 0, 2}, + { 2, 0, 1}, + { 2, 1, 0}}}; - for (int pb = 0; pb < 2; ++pb){ - for (int pg = 0; pg < 2; ++pg){ - for (int qa = 0; qa < 2; ++qa){ - for (int qg = 0; qg < 2; ++qg){ - for (int ra = 0; ra < 2; ++ra){ - for (int rb = 0; rb < 2; ++rb){ - const double cprod = C(vect[0], 1, pb, pg)*C(vect[1], qa, 1, qg)*C(vect[2], ra, rb, 1); - const double denom = (qa+ra+1) * (pb+rb+1) * (pg+qg+1); - volume += perm123sign(permutation[0], permutation[1], permutation[2]) * cprod / denom; - } - } - } - } - } - } - } while (std::next_permutation(std::begin(permutation), std::end(permutation))); + static const std::array pqr_array + = {{{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 1, 1}, + {0, 0, 0, 1, 0, 0}, {0, 0, 0, 1, 0, 1}, {0, 0, 0, 1, 1, 0}, {0, 0, 0, 1, 1, 1}, + {0, 0, 1, 0, 0, 0}, {0, 0, 1, 0, 0, 1}, {0, 0, 1, 0, 1, 0}, {0, 0, 1, 0, 1, 1}, + {0, 0, 1, 1, 0, 0}, {0, 0, 1, 1, 0, 1}, {0, 0, 1, 1, 1, 0}, {0, 0, 1, 1, 1, 1}, + {0, 1, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 1}, {0, 1, 0, 0, 1, 0}, {0, 1, 0, 0, 1, 1}, + {0, 1, 0, 1, 0, 0}, {0, 1, 0, 1, 0, 1}, {0, 1, 0, 1, 1, 0}, {0, 1, 0, 1, 1, 1}, + {0, 1, 1, 0, 0, 0}, {0, 1, 1, 0, 0, 1}, {0, 1, 1, 0, 1, 0}, {0, 1, 1, 0, 1, 1}, + {0, 1, 1, 1, 0, 0}, {0, 1, 1, 1, 0, 1}, {0, 1, 1, 1, 1, 0}, {0, 1, 1, 1, 1, 1}, + {1, 0, 0, 0, 0, 0}, {1, 0, 0, 0, 0, 1}, {1, 0, 0, 0, 1, 0}, {1, 0, 0, 0, 1, 1}, + {1, 0, 0, 1, 0, 0}, {1, 0, 0, 1, 0, 1}, {1, 0, 0, 1, 1, 0}, {1, 0, 0, 1, 1, 1}, + {1, 0, 1, 0, 0, 0}, {1, 0, 1, 0, 0, 1}, {1, 0, 1, 0, 1, 0}, {1, 0, 1, 0, 1, 1}, + {1, 0, 1, 1, 0, 0}, {1, 0, 1, 1, 0, 1}, {1, 0, 1, 1, 1, 0}, {1, 0, 1, 1, 1, 1}, + {1, 1, 0, 0, 0, 0}, {1, 1, 0, 0, 0, 1}, {1, 1, 0, 0, 1, 0}, {1, 1, 0, 0, 1, 1}, + {1, 1, 0, 1, 0, 0}, {1, 1, 0, 1, 0, 1}, {1, 1, 0, 1, 1, 0}, {1, 1, 0, 1, 1, 1}, + {1, 1, 1, 0, 0, 0}, {1, 1, 1, 0, 0, 1}, {1, 1, 1, 0, 1, 0}, {1, 1, 1, 0, 1, 1}, + {1, 1, 1, 1, 0, 0}, {1, 1, 1, 1, 0, 1}, {1, 1, 1, 1, 1, 0}, {1, 1, 1, 1, 1, 1}}}; - return std::fabs(volume); + double volume = 0.0; + const double* vect[3]; + const std::array,3> data = {{X, Y, Z}}; + double perm_sign = 1; + for (const auto& perm : permutation) { + for (std::size_t perm_index = 0; perm_index < 3; perm_index++) + vect[perm_index] = data[perm[perm_index]].data(); + + for (const auto& pqr : pqr_array) { + const double cprod = C(vect[0], 1, pqr.pb, pqr.pg)*C(vect[1], pqr.qa, 1, pqr.qg)*C(vect[2], pqr.ra, pqr.rb, 1); + const double denom = (pqr.qa + pqr.ra + 1) * (pqr.pb + pqr.rb + 1) * (pqr.pg + pqr.qg + 1); + volume += perm_sign * cprod / denom; + } + + perm_sign *= -1; + } + return std::fabs(volume); }