From 27e88a54f268dae79529b16e6014b4d8862f14e8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Fri, 2 Jan 2015 15:22:02 +0100
Subject: [PATCH] Fix anisotropy closeness test.

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.
---
 opm/core/tof/AnisotropicEikonal.cpp | 33 +++++++++++++++++++++++------
 opm/core/tof/AnisotropicEikonal.hpp |  7 +++++-
 2 files changed, 33 insertions(+), 7 deletions(-)

diff --git a/opm/core/tof/AnisotropicEikonal.cpp b/opm/core/tof/AnisotropicEikonal.cpp
index 16c3f0fe..89e04461 100644
--- a/opm/core/tof/AnisotropicEikonal.cpp
+++ b/opm/core/tof/AnisotropicEikonal.cpp
@@ -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
diff --git a/opm/core/tof/AnisotropicEikonal.hpp b/opm/core/tof/AnisotropicEikonal.hpp
index 94daedb9..942c4d97 100644
--- a/opm/core/tof/AnisotropicEikonal.hpp
+++ b/opm/core/tof/AnisotropicEikonal.hpp
@@ -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