From 0f861e64ed1a31f318a8d813ea86c596406ec2ca Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Thu, 20 Dec 2012 18:03:08 +0100
Subject: [PATCH] Limiter now works reasonably well.

Case with no inflow faces should be checked.
---
 .../TransportModelTracerTofDiscGal.cpp        | 41 ++++++++++---------
 1 file changed, 22 insertions(+), 19 deletions(-)

diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
index 9738dc9c..0dad656e 100644
--- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
+++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
@@ -404,9 +404,10 @@ namespace Opm
         const int dim = grid_.dimensions;
         const int num_basis = DGBasis::numBasisFunc(dim, degree_);
 
-        double limiter = 1e100;
+        double max_slope_mult = 0.0;
+        int num_upstream_faces = 0;
         // For inflow faces, ensure that cell tof does not dip below
-        // the minimum value from upstream (for that face).
+        // the minimum value from upstream (for all faces).
         for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
             const int face = grid_.cell_faces[hface];
             double flux = 0.0;
@@ -418,6 +419,11 @@ namespace Opm
                 flux = -darcyflux_[face];
                 upstream_cell = grid_.face_cells[2*face];
             }
+            if (flux >= 0.0) {
+                // This is a downstream face.
+                continue;
+            }
+            ++num_upstream_faces;
 
             // Evaluate the solution in all corners, and find the appropriate limiter.
             bool upstream = (upstream_cell >= 0 && flux < 0.0);
@@ -437,26 +443,23 @@ namespace Opm
                     min_upstream = std::min(min_upstream, tof_upstream);
                 }
             }
-            if (min_here < min_upstream) {
-                // Must limit slope.
-                const double tof_c = tof_coeff_[num_basis*cell];
-                if (tof_c < min_upstream) {
-                    // Handle by setting a flat solution.
-                    std::cout << "Trouble in cell " << cell << std::endl;
-                    limiter = 0.0;
-                    tof_coeff_[num_basis*cell] = min_upstream;
-                    break;
-                }
-                const double face_limit = (tof_c - min_upstream)/(tof_c - min_here);
-                limiter = std::min(limiter, face_limit);
+            // Compute maximum slope multiplier.
+            const double tof_c = tof_coeff_[num_basis*cell];
+            if (tof_c < min_upstream) {
+                // Handle by setting a flat solution.
+                std::cout << "Trouble in cell " << cell << std::endl;
+                max_slope_mult = 0.0;
+                tof_coeff_[num_basis*cell] = min_upstream;
+                break;
             }
+            const double face_mult = (tof_c - min_upstream)/(tof_c - min_here);
+            max_slope_mult = std::max(max_slope_mult, face_mult);
         }
+        ASSERT(max_slope_mult >= 0.0);
 
-        if (limiter < 0.0) {
-            THROW("Error in limiter.");
-        }
-
-        if (limiter < 1.0) {
+        // Actually do the limiting (if applicable).
+        const double limiter = max_slope_mult;
+        if (num_upstream_faces > 0 && limiter < 1.0) {
             std::cout << "Applying limiter in cell " << cell << ", limiter = " << limiter << std::endl;
             for (int i = num_basis*cell + 1; i < num_basis*(cell+1); ++i) {
                 tof_coeff_[i] *= limiter;