diff --git a/opm/core/tof/AnisotropicEikonal.cpp b/opm/core/tof/AnisotropicEikonal.cpp index 2ce60637b..e99dc888a 100644 --- a/opm/core/tof/AnisotropicEikonal.cpp +++ b/opm/core/tof/AnisotropicEikonal.cpp @@ -20,7 +20,7 @@ #include #include #include -#include +#include namespace Opm { @@ -92,7 +92,7 @@ namespace Opm for (int nb = 0; nb < num_nb; ++nb) { const int nb_cell = cell_neighbours_[scell][nb]; if (!is_accepted_[nb_cell] && !is_considered_[nb_cell]) { - const double value = computeValue(nb_cell, metric); + const double value = computeValue(nb_cell, metric, solution.data()); pushConsidered(std::make_pair(value, nb_cell)); } } @@ -130,7 +130,7 @@ namespace Opm const int nb_cell = *it; if (!is_accepted_[nb_cell] && !is_considered_[nb_cell]) { assert(solution[nb_cell] == inf); - const double value = computeValue(nb_cell, metric); + const double value = computeValue(nb_cell, metric, solution.data()); pushConsidered(std::make_pair(value, nb_cell)); } } @@ -140,7 +140,7 @@ namespace Opm for (auto it = considered_.begin(); it != considered_.end(); ++it) { const int ccell = it->second; if (isClose(rcell, ccell, metric)) { - const double value = computeValue(ccell, metric); + const double value = computeValue(ccell, metric, solution.data()); if (value < it->first) { // Update value for considered cell. // Note that as solution values decrease, their @@ -173,7 +173,8 @@ namespace Opm double AnisotropicEikonal2d::computeValue(const int cell, - const double* metric) const + const double* metric, + const double* solution) const { const auto& nbs = cell_neighbours_[cell]; const int num_nbs = nbs.size(); @@ -182,7 +183,7 @@ namespace Opm for (int ii = 0; ii < num_nbs; ++ii) { const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] }; if (accepted_front_.count(n[0]) && accepted_front_.count(n[1])) { - const double cand_val = computeFromTri(cell, n[0], n[1], metric); + const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution); val = std::min(val, cand_val); } } @@ -191,7 +192,7 @@ namespace Opm // so we go for a single-neighbour update. for (int ii = 0; ii < num_nbs; ++ii) { if (accepted_front_.count(nbs[ii])) { - const double cand_val = computeFromLine(cell, nbs[ii], metric); + const double cand_val = computeFromLine(cell, nbs[ii], metric, solution); val = std::min(val, cand_val); } } @@ -204,28 +205,85 @@ namespace Opm + double distanceAniso(const double v1[2], + const double v2[2], + const double g[4]) + { + const double d[2] = { v2[0] - v1[0], v2[1] - v1[1] }; + const double dist = std::sqrt(+ g[0] * d[0] * d[0] + + g[1] * d[0] * d[1] + + g[2] * d[1] * d[0] + + g[3] * d[1] * d[1]); + return dist; + } + + + + + double AnisotropicEikonal2d::computeFromLine(const int cell, const int from, - const double* metric) const + const double* metric, + const double* solution) const { assert(!is_accepted_[cell]); assert(is_accepted_[from]); - return 2.0; + // Applying the first fundamental form to compute geodesic distance. + // Using the metric of 'cell', not 'from'. + const double dist = distanceAniso(grid_.cell_centroids + 2 * cell, + grid_.cell_centroids + 2 * from, + metric + 4 * cell); + return solution[from] + dist; } + struct DistanceDerivative + { + const double* x1; + const double* x2; + const double* x; + double u1; + double u2; + const double* g; + double operator()(const double theta) const + { + double a[2] = { x[0] - (1-theta)*x1[0] + theta*x2[0], x[1] - (1-theta)*x1[1] + theta*x2[1] }; + double b[2] = { x1[0] - x2[0], x1[1] - x2[1] }; + return u2 - u1 + 2*(a[0]*b[0]*g[0] + a[0]*b[1]*g[1] + a[1]*b[0]*g[2] + a[1]*b[1]*g[3]); + } + }; + + + + + double AnisotropicEikonal2d::computeFromTri(const int cell, const int n0, const int n1, - const double* metric) const + const double* metric, + const double* solution) const { assert(!is_accepted_[cell]); assert(is_accepted_[n0]); assert(is_accepted_[n1]); - return 1.0; + DistanceDerivative dd; + dd.x1 = grid_.cell_centroids + 2 * n0; + dd.x2 = grid_.cell_centroids + 2 * n1; + dd.x = grid_.cell_centroids + 2 * cell; + dd.u1 = solution[n0]; + dd.u2 = solution[n1]; + dd.g = metric + 4 * cell; + int iter = 0; + const double theta = RegulaFalsi::solve(dd, 0.0, 1.0, 15, 1e-8, iter); + const double xt[2] = { (1-theta)*dd.x1[0] + theta*dd.x2[0], + (1-theta)*dd.x1[1] + theta*dd.x2[1] }; + const double d1 = distanceAniso(dd.x1, dd.x, dd.g) + solution[n0]; + const double d2 = distanceAniso(dd.x2, dd.x, dd.g) + solution[n1]; + const double dt = distanceAniso(xt, dd.x, dd.g) + (1-theta)*solution[n0] + theta*solution[n1]; + return std::min(d1, std::min(d2, dt)); } diff --git a/opm/core/tof/AnisotropicEikonal.hpp b/opm/core/tof/AnisotropicEikonal.hpp index b72a8b8ef..c977b0db9 100644 --- a/opm/core/tof/AnisotropicEikonal.hpp +++ b/opm/core/tof/AnisotropicEikonal.hpp @@ -66,9 +66,9 @@ namespace Opm std::vector is_considered_; bool isClose(const int c1, const int c2, const double* metric) const; - double computeValue(const int cell, const double* metric) const; - double computeFromLine(const int cell, const int from, const double* metric) const; - double computeFromTri(const int cell, const int n0, const int n1, const double* metric) const; + double computeValue(const int cell, const double* metric, const double* solution) const; + double computeFromLine(const int cell, const int from, const double* metric, const double* solution) const; + double computeFromTri(const int cell, const int n0, const int n1, const double* metric, const double* solution) const; const ValueAndCell& topConsidered() const; void pushConsidered(const ValueAndCell& vc);