OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-27 01:31:00 -06:00
Code Issues Packages Projects Releases Wiki Activity

Whitespace fix.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2014-12-01 15:17:30 +01:00
parent 451bca1753
commit 2eab160761
2 changed files with 88 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

144
opm/core/tof/AnisotropicEikonal.cpp
View File

@ -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<int>& startcells,
std::vector<double>& solution)
const std::vector<int>& startcells,
std::vector<double>& 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();
}

32
opm/core/tof/AnisotropicEikonal.hpp
View File

@ -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<int>& startcells,
std::vector<double>& solution);
const std::vector<int>& startcells,
std::vector<double>& solution);
private:
// Grid and topology.
const UnstructuredGrid& grid_;
SparseTable<int> cell_neighbours_;
const UnstructuredGrid& grid_;
SparseTable<int> cell_neighbours_;
// Keep track of accepted cells.
std::vector<char> is_accepted_;
std::vector<char> is_accepted_;
std::set<int> accepted_front_;
// Keep track of considered cells.
typedef std::pair<double, int> ValueAndCell;
typedef std::pair<double, int> ValueAndCell;
typedef boost::heap::compare<std::greater<ValueAndCell>> Comparator;
typedef boost::heap::fibonacci_heap<ValueAndCell, Comparator> Heap;
Heap considered_;
typedef Heap::handle_type HeapHandle;
std::map<int, HeapHandle> considered_handles_;
std::vector<char> is_considered_;
std::vector<char> 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