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

Fix anisotropy closeness test.

Browse Source 
The old test was simply wrong: it computed the M-distance and compared
to the grid radius, which becomes dependent on the scaling of the
metric M. The corrected test in isClose() depends on the anisotropy
ratio of M and the grid radius.
This commit is contained in:
Atgeirr Flø Rasmussen 2015-01-02 15:22:02 +01:00
parent 8491d186d4
commit 1e39010ba4
2 changed files with 33 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -59,7 +59,8 @@ namespace Opm
/// Construct solver.
/// \param[in] grid A 2d grid.
AnisotropicEikonal2d::AnisotropicEikonal2d(const UnstructuredGrid& grid)
: grid_(grid)
: grid_(grid),
safety_factor_(1.2)
{
if (grid.dimensions != 2) {
OPM_THROW(std::logic_error, "Grid for AnisotropicEikonal2d must be 2d.");
@ -77,6 +78,9 @@ namespace Opm
const std::vector<int>& startcells,
std::vector<double>& solution)
{
// Compute anisotropy ratios to be used by isClose().
computeAnisoRatio(metric);
// 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,
@ -162,7 +166,7 @@ namespace Opm
// distance h * F_2/F1 from x_r. Use min of previous and new.
for (auto it = considered_.begin(); it != considered_.end(); ++it) {
const int ccell = it->second;
if (isClose(rcell, ccell, metric)) {
if (isClose(rcell, ccell)) {
const double value = computeValueUpdate(ccell, metric, solution.data(), rcell);
if (value < it->first) {
// Update value for considered cell.
@ -195,13 +199,11 @@ namespace Opm
bool AnisotropicEikonal2d::isClose(const int c1,
const int c2,
const double* metric) const
const int c2) const
{
const double* v[] = { grid_.cell_centroids + 2*c1,
grid_.cell_centroids + 2*c2 };
const double* m = metric + 4*c1;
return distanceAniso(v[0], v[1], m) < 3.0 * grid_radius_[c1];
return distanceIso(v[0], v[1]) < safety_factor_ * aniso_ratio_[c1] * grid_radius_[c1];
}
@ -402,5 +404,24 @@ namespace Opm
void AnisotropicEikonal2d::computeAnisoRatio(const double* metric)
{
const int num_cells = cell_neighbours_.size();
aniso_ratio_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
const double* m = metric + 4*cell;
// Find the two eigenvalues from trace and determinant.
const double t = m[0] + m[3];
const double d = m[0]*m[3] - m[1]*m[2];
const double sd = std::sqrt(t*t/4.0 - d);
const double eig[2] = { t/2.0 - sd, t/2.0 + sd };
// Anisotropy ratio is the max ratio of the eigenvalues.
aniso_ratio_[cell] = std::max(eig[0]/eig[1], eig[1]/eig[0]);
}
}
} // namespace Opm

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

@ -55,7 +55,11 @@ namespace Opm
// Keep track of accepted cells.
std::vector<char> is_accepted_;
std::set<int> accepted_front_;
// Quantities relating to anisotropy.
std::vector<double> grid_radius_;
std::vector<double> aniso_ratio_;
const double safety_factor_;
// Keep track of considered cells.
typedef std::pair<double, int> ValueAndCell;
@ -66,7 +70,7 @@ namespace Opm
std::map<int, HeapHandle> considered_handles_;
std::vector<char> is_considered_;
bool isClose(const int c1, const int c2, const double* metric) const;
bool isClose(const int c1, const int c2) 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;
@ -77,6 +81,7 @@ namespace Opm
void popConsidered();
void computeGridRadius();
void computeAnisoRatio(const double* metric);
};
} // namespace Opm