From 15ab570cd5b38595c6534d8dab9b331f9d27147b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Wed, 24 Apr 2013 11:27:04 +0200
Subject: [PATCH] Added limiter for tracer.

---
 opm/core/tof/TofDiscGalReorder.cpp | 48 +++++++++++++++++++++++++++++-
 opm/core/tof/TofDiscGalReorder.hpp |  4 +++
 2 files changed, 51 insertions(+), 1 deletion(-)

diff --git a/opm/core/tof/TofDiscGalReorder.cpp b/opm/core/tof/TofDiscGalReorder.cpp
index e2e4839a6..95ef468bf 100644
--- a/opm/core/tof/TofDiscGalReorder.cpp
+++ b/opm/core/tof/TofDiscGalReorder.cpp
@@ -500,7 +500,11 @@ namespace Opm
             std::cout << std::endl;
 #endif
             applyLimiter(cell, tof_coeff_);
-            // We do not (yet) apply a limiter to the tracer solution.
+            if (num_tracers_ && tracerhead_by_cell_[cell] == NoTracerHead) {
+                for (int tr = 0; tr < num_tracers_; ++tr) {
+                    applyTracerLimiter(cell, tracer_coeff_ + cell*num_tracers_*num_basis + tr*num_basis);
+                }
+            }
         }
     }
 
@@ -721,5 +725,47 @@ namespace Opm
 
 
 
+    void TofDiscGalReorder::applyTracerLimiter(const int cell, double* local_coeff)
+    {
+        // Evaluate the solution in all corners of all faces. Extract max and min.
+        const int dim = grid_.dimensions;
+        const int num_basis = basis_func_->numBasisFunc();
+        double min_cornerval = 1e100;
+        double max_cornerval = -1e100;
+        for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
+            const int face = grid_.cell_faces[hface];
+            for (int fnode = grid_.face_nodepos[face]; fnode < grid_.face_nodepos[face+1]; ++fnode) {
+                const double* nc = grid_.node_coordinates + dim*grid_.face_nodes[fnode];
+                basis_func_->eval(cell, nc, &basis_[0]);
+                const double tracer_corner = std::inner_product(basis_.begin(), basis_.end(),
+                                                                local_coeff, 0.0);
+                min_cornerval = std::min(min_cornerval, tracer_corner);
+                max_cornerval = std::max(min_cornerval, tracer_corner);
+            }
+        }
+        const double average = basis_func_->functionAverage(local_coeff);
+        if (average < 0.0 || average > 1.0) {
+            // Adjust average. Flatten gradient.
+            std::fill(local_coeff, local_coeff + num_basis, 0.0);
+            if (average > 1.0) {
+                basis_func_->addConstant(1.0, local_coeff);
+            }
+        } else {
+            // Possibly adjust gradient.
+            double factor = 1.0;
+            if (min_cornerval < 0.0) {
+                factor = average/(average - min_cornerval);
+            }
+            if (max_cornerval > 1.0) {
+                factor = std::min(factor, (1.0 - average)/(max_cornerval - average));
+            }
+            if (factor != 1.0) {
+                basis_func_->multiplyGradient(factor, local_coeff);
+            }
+        }
+    }
+
+
+
 
 } // namespace Opm
diff --git a/opm/core/tof/TofDiscGalReorder.hpp b/opm/core/tof/TofDiscGalReorder.hpp
index c4e67f2df..b13affac8 100644
--- a/opm/core/tof/TofDiscGalReorder.hpp
+++ b/opm/core/tof/TofDiscGalReorder.hpp
@@ -174,6 +174,10 @@ namespace Opm
         void applyLimiterAsSimultaneousPostProcess();
         double totalFlux(const int cell) const;
         double minCornerVal(const int cell, const int face) const;
+
+        // Apply a simple (restrict to [0,1]) limiter.
+        // Intended for tracers.
+        void applyTracerLimiter(const int cell, double* local_coeff);
     };
 
 } // namespace Opm