Added solveMultiCell() interface and implementation.

- Added solveMultiCell() virtual method.
  - TransportModelInterface::reorderAndTransport() now calls solveMultiCell()
    instead of aborting if encountering multi-cell components.
  - Implemented solveMultiCell() in TransportModelTwophase by solving
    each cell individually with solveSingleCell() and repeating until
    saturation change is small (hardcoded 1e-9 for now).

opm/core/transport/reorder/TransportModelInterface.cpp

@@ -27,22 +27,27 @@
 
 void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux)
 {
-    // Compute sequence of single-cell problems
+    // Compute reordered sequence of single-cell problems
     std::vector<int> sequence(grid.number_of_cells);
     std::vector<int> components(grid.number_of_cells + 1);
     int ncomponents;
     compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents);
 
-    // Assume all strong components are single-cell domains.
+    // Invoke appropriate solve method for each interdependent component.
-    assert(ncomponents == grid.number_of_cells);
+    for (int comp = 0; comp < ncomponents; ++comp) {
-    for (int i = 0; i < grid.number_of_cells; ++i) {
 #ifdef MATLAB_MEX_FILE
+	// \TODO replace this with general signal handling code, check if it costs performance.
         if (interrupt_signal) {
             mexPrintf("Reorder loop interrupted by user: %d of %d "
                       "cells finished.\n", i, grid.number_of_cells);
             break;
         }
 #endif
-        solveSingleCell(sequence[i]);
+	const int comp_size = components[comp + 1] - components[comp];
+	if (comp_size == 1) {
+	    solveSingleCell(sequence[components[comp]]);
+	} else {
+	    solveMultiCell(comp_size, &sequence[components[comp]]);
+	}
     }
 }

opm/core/transport/reorder/TransportModelInterface.hpp

@@ -26,17 +26,18 @@ namespace Opm
 {
 
     /// Interface for reordering transport models.
-    /// A transport model must provide the solveSingleCell()
+    /// A transport model must provide the solveSingleCell() and
-    /// method, and is expected to implement a solve() method
+    /// solveMultiCell methods, and is expected to implement a solve()
-    /// that will have an interface geared to the model's needs.
+    /// method that will have an interface geared to the model's
-    /// (The solve() method is therefore not virtual in this class).
+    /// needs. (The solve() method is therefore not virtual in this
-    /// The reorderAndTransport() method is provided as an
+    /// class.) The reorderAndTransport() method is provided as an aid
-    /// aid to implementing solve() in subclasses.
+    /// to implementing solve() in subclasses.
     class TransportModelInterface
     {
     public:
 	virtual ~TransportModelInterface() {}
-	virtual void solveSingleCell(int cell) = 0;
+	virtual void solveSingleCell(const int cell) = 0;
+	virtual void solveMultiCell(const int num_cells, const int* cells) = 0;
     protected:
 	void reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux);
     };

opm/core/transport/reorder/TransportModelTwophase.cpp

@@ -23,6 +23,9 @@
 #include <opm/core/transport/reorder/nlsolvers.h>
 #include <opm/core/utility/RootFinders.hpp>
 
+#include <fstream>
+#include <iterator>
+
 namespace Opm
 {
 
@@ -111,18 +114,64 @@ namespace Opm
     };
 
 
-    void TransportModelTwophase::solveSingleCell(int cell)
+    void TransportModelTwophase::solveSingleCell(const int cell)
     {
 	Residual res(*this, cell);
+	const double tol = 1e-9;
+	// const double r0 = res(saturation_[cell]);
+	// if (std::fabs(r0) < tol) {
+	//     return;
+	// }
 	const double a = 0.0;
 	const double b = 1.0;
 	const int maxit = 20;
-	const double tol = 1e-9;
 	int iters_used;
+	// saturation_[cell] = modifiedRegulaFalsi(res, a, b, saturation_[cell], maxit, tol, iters_used);
 	saturation_[cell] = modifiedRegulaFalsi(res, a, b, maxit, tol, iters_used);
 	fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
     }
 
+    void TransportModelTwophase::solveMultiCell(const int num_cells, const int* cells)
+    {
+	// std::ofstream os("dump");
+	// std::copy(cells, cells + num_cells, std::ostream_iterator<double>(os, "\n"));
+	double max_s_change = 0.0;
+	const double tol = 1e-9;
+	const int max_iters = 300;
+	int num_iters = 0;
+	// Must store s0 before we start.
+	std::vector<double> s0(num_cells);
+	// Must set initial fractional flows before we start.
+	for (int i = 0; i < num_cells; ++i) {
+	    const int cell = cells[i];
+	    fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
+	    s0[i] = saturation_[cell];
+	}
+	do {
+	    int max_change_cell = -1;
+	    max_s_change = 0.0;
+	    for (int i = 0; i < num_cells; ++i) {
+		const int cell = cells[i];
+		const double old_s = saturation_[cell];
+		saturation_[cell] = s0[i];
+		solveSingleCell(cell);
+		// std::cout << "delta s = " << saturation_[cell] - old_s << std::endl;
+		if (max_s_change < std::fabs(saturation_[cell] - old_s)) {
+		    max_change_cell = cell;
+		}
+		max_s_change = std::max(max_s_change, std::fabs(saturation_[cell] - old_s));
+	    }
+	    // std::cout << "Iter = " << num_iters << "    max_s_change = " << max_s_change
+	    // 	      << "    in cell " << max_change_cell << std::endl;
+	} while (max_s_change > tol && ++num_iters < max_iters);
+	if (max_s_change > tol) {
+	    THROW("In solveMultiCell(), we did not converge after "
+	    	  << num_iters << " iterations. Delta s = " << max_s_change);
+	}
+	std::cout << "Solved " << num_cells << " cell multicell problem in "
+		  << num_iters << " iterations." << std::endl;
+    }
+
     double TransportModelTwophase::fracFlow(double s, int cell) const
     {
 	double sat[2] = { s, 1.0 - s };

opm/core/transport/reorder/TransportModelTwophase.hpp

@@ -42,7 +42,8 @@ namespace Opm
 		   const double dt,
 		   double* saturation);
 
-	virtual void solveSingleCell(int cell);
+	virtual void solveSingleCell(const int cell);
+	virtual void solveMultiCell(const int num_cells, const int* cells);
 
     private:
 	const UnstructuredGrid& grid_;