diff --git a/opm/core/tof/AnisotropicEikonal.cpp b/opm/core/tof/AnisotropicEikonal.cpp index 0840aaa6..f715e36a 100644 --- a/opm/core/tof/AnisotropicEikonal.cpp +++ b/opm/core/tof/AnisotropicEikonal.cpp @@ -28,13 +28,13 @@ namespace Opm /// Construct solver. /// \param[in] grid A 2d grid. AnisotropicEikonal2d::AnisotropicEikonal2d(const UnstructuredGrid& grid) - : grid_(grid) + : grid_(grid) { - if (grid.dimensions != 2) { - OPM_THROW(std::logic_error, "Grid for AnisotropicEikonal2d must be 2d."); - } - cell_neighbours_ = cellNeighboursAcrossVertices(grid); - orderCounterClockwise(grid, cell_neighbours_); + if (grid.dimensions != 2) { + OPM_THROW(std::logic_error, "Grid for AnisotropicEikonal2d must be 2d."); + } + cell_neighbours_ = cellNeighboursAcrossVertices(grid); + orderCounterClockwise(grid, cell_neighbours_); } /// Solve the eikonal equation. @@ -42,61 +42,61 @@ namespace Opm /// \param[in] startcells Array of cells where u = 0 at the centroid. /// \param[out] solution Array of solution to the eikonal equation. void AnisotropicEikonal2d::solve(const double* metric, - const std::vector& startcells, - std::vector& solution) + const std::vector& startcells, + std::vector& solution) { - // The algorithm used is described in J.A. Sethian and A. Vladimirsky, - // "Ordered Upwind Methods for Static Hamilton-Jacobi Equations". - // Notation in comments is as used in that paper: U is the solution, - // and q is the boundary condition. One difference is that we talk about - // grid cells instead of mesh points. - // - // Algorithm summary: - // 1. Put all cells in Far. U_i = \inf. - // 2. Move the startcells to Accepted. U_i = q(x_i) - // 3. Move cells adjacent to startcells to Considered, evaluate - // U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k} - // 4. Find the Considered cell with the smallest value: r. - // 5. Move cell r to Accepted. Update AcceptedFront. - // 6. Recompute the value for all Considered cells within - // distance h * F_2/F1 from x_r. Use min of previous and new. - // 7. Move cells adjacent to r from Far to Considered. - // 8. If Considered is not empty, go to step 4. + // The algorithm used is described in J.A. Sethian and A. Vladimirsky, + // "Ordered Upwind Methods for Static Hamilton-Jacobi Equations". + // Notation in comments is as used in that paper: U is the solution, + // and q is the boundary condition. One difference is that we talk about + // grid cells instead of mesh points. + // + // Algorithm summary: + // 1. Put all cells in Far. U_i = \inf. + // 2. Move the startcells to Accepted. U_i = q(x_i) + // 3. Move cells adjacent to startcells to Considered, evaluate + // U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k} + // 4. Find the Considered cell with the smallest value: r. + // 5. Move cell r to Accepted. Update AcceptedFront. + // 6. Recompute the value for all Considered cells within + // distance h * F_2/F1 from x_r. Use min of previous and new. + // 7. Move cells adjacent to r from Far to Considered. + // 8. If Considered is not empty, go to step 4. - // 1. Put all cells in Far. U_i = \inf. - const int num_cells = grid_.number_of_cells; - const double inf = 1e100; - solution.clear(); - solution.resize(num_cells, inf); - is_accepted_.clear(); - is_accepted_.resize(num_cells, false); + // 1. Put all cells in Far. U_i = \inf. + const int num_cells = grid_.number_of_cells; + const double inf = 1e100; + solution.clear(); + solution.resize(num_cells, inf); + is_accepted_.clear(); + is_accepted_.resize(num_cells, false); accepted_front_.clear(); - considered_.clear(); + considered_.clear(); considered_handles_.clear(); - is_considered_.clear(); - is_considered_.resize(num_cells, false); + is_considered_.clear(); + is_considered_.resize(num_cells, false); - // 2. Move the startcells to Accepted. U_i = q(x_i) - const int num_startcells = startcells.size(); - for (int ii = 0; ii < num_startcells; ++ii) { - is_accepted_[startcells[ii]] = true; - solution[startcells[ii]] = 0.0; - } - accepted_front_.insert(startcells.begin(), startcells.end()); + // 2. Move the startcells to Accepted. U_i = q(x_i) + const int num_startcells = startcells.size(); + for (int ii = 0; ii < num_startcells; ++ii) { + is_accepted_[startcells[ii]] = true; + solution[startcells[ii]] = 0.0; + } + accepted_front_.insert(startcells.begin(), startcells.end()); - // 3. Move cells adjacent to startcells to Considered, evaluate - // U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k} - for (int ii = 0; ii < num_startcells; ++ii) { - const int scell = startcells[ii]; - const int num_nb = cell_neighbours_[scell].size(); - 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, solution.data()); - pushConsidered(std::make_pair(value, nb_cell)); - } - } - } + // 3. Move cells adjacent to startcells to Considered, evaluate + // U_i = min_{(x_j,x_k) \in NF(x_i)} G_{j,k} + for (int ii = 0; ii < num_startcells; ++ii) { + const int scell = startcells[ii]; + const int num_nb = cell_neighbours_[scell].size(); + 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, solution.data()); + pushConsidered(std::make_pair(value, nb_cell)); + } + } + } while (!considered_.empty()) { // 4. Find the Considered cell with the smallest value: r. @@ -178,17 +178,17 @@ namespace Opm const double* solution) const { // std::cout << "++++ computeValue(), cell = " << cell << std::endl; - const auto& nbs = cell_neighbours_[cell]; - const int num_nbs = nbs.size(); + const auto& nbs = cell_neighbours_[cell]; + const int num_nbs = nbs.size(); const double inf = 1e100; - double val = inf; - for (int ii = 0; ii < num_nbs; ++ii) { - const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] }; + double val = inf; + 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, solution); val = std::min(val, cand_val); } - } + } if (val == inf) { // Failed to find two accepted front nodes adjacent to this, // so we go for a single-neighbour update. @@ -201,7 +201,7 @@ namespace Opm } assert(val != inf); // std::cout << "---> " << val << std::endl; - return val; + return val; } @@ -214,18 +214,18 @@ namespace Opm const int new_cell) const { // std::cout << "++++ computeValueUpdate(), cell = " << cell << std::endl; - const auto& nbs = cell_neighbours_[cell]; - const int num_nbs = nbs.size(); + const auto& nbs = cell_neighbours_[cell]; + const int num_nbs = nbs.size(); const double inf = 1e100; - double val = inf; - for (int ii = 0; ii < num_nbs; ++ii) { - const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] }; + double val = inf; + for (int ii = 0; ii < num_nbs; ++ii) { + const int n[2] = { nbs[ii], nbs[(ii+1) % num_nbs] }; if ((n[0] == new_cell || n[1] == new_cell) && accepted_front_.count(n[0]) && accepted_front_.count(n[1])) { const double cand_val = computeFromTri(cell, n[0], n[1], metric, solution); val = std::min(val, cand_val); } - } + } if (val == inf) { // Failed to find two accepted front nodes adjacent to this, // so we go for a single-neighbour update. @@ -237,7 +237,7 @@ namespace Opm } } // std::cout << "---> " << val << std::endl; - return val; + return val; } @@ -336,7 +336,7 @@ namespace Opm const AnisotropicEikonal2d::ValueAndCell& AnisotropicEikonal2d::topConsidered() const { - return considered_.top(); + return considered_.top(); } @@ -345,7 +345,7 @@ namespace Opm void AnisotropicEikonal2d::pushConsidered(const ValueAndCell& vc) { - HeapHandle h = considered_.push(vc); + HeapHandle h = considered_.push(vc); considered_handles_[vc.second] = h; is_considered_[vc.second] = true; } @@ -358,7 +358,7 @@ namespace Opm { is_considered_[considered_.top().second] = false; considered_handles_.erase(considered_.top().second); - considered_.pop(); + considered_.pop(); } diff --git a/opm/core/tof/AnisotropicEikonal.hpp b/opm/core/tof/AnisotropicEikonal.hpp index 1e2480c2..215219e4 100644 --- a/opm/core/tof/AnisotropicEikonal.hpp +++ b/opm/core/tof/AnisotropicEikonal.hpp @@ -36,44 +36,44 @@ namespace Opm class AnisotropicEikonal2d { public: - /// Construct solver. + /// Construct solver. /// \param[in] grid A 2d grid. explicit AnisotropicEikonal2d(const UnstructuredGrid& grid); /// Solve the eikonal equation. - /// \param[in] metric Array of metric tensors, M, for each cell. + /// \param[in] metric Array of metric tensors, M, for each cell. /// \param[in] startcells Array of cells where u = 0 at the centroid. /// \param[out] solution Array of solution to the eikonal equation. void solve(const double* metric, - const std::vector& startcells, - std::vector& solution); + const std::vector& startcells, + std::vector& solution); private: // Grid and topology. - const UnstructuredGrid& grid_; - SparseTable cell_neighbours_; + const UnstructuredGrid& grid_; + SparseTable cell_neighbours_; // Keep track of accepted cells. - std::vector is_accepted_; + std::vector is_accepted_; std::set accepted_front_; // Keep track of considered cells. - typedef std::pair ValueAndCell; + typedef std::pair ValueAndCell; typedef boost::heap::compare> Comparator; typedef boost::heap::fibonacci_heap Heap; Heap considered_; typedef Heap::handle_type HeapHandle; std::map considered_handles_; - std::vector is_considered_; + 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* solution) const; - double computeValueUpdate(const int cell, const double* metric, const double* solution, const int new_cell) 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; + double computeValue(const int cell, const double* metric, const double* solution) const; + double computeValueUpdate(const int cell, const double* metric, const double* solution, const int new_cell) 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); - void popConsidered(); + const ValueAndCell& topConsidered() const; + void pushConsidered(const ValueAndCell& vc); + void popConsidered(); }; } // namespace Opm