Merge from upstream.

examples/spu_2p.cpp

@@ -368,7 +368,7 @@ main(int argc, char** argv)
     if (use_column_solver) {
         Opm::extractColumn(*grid->c_grid(), columns);
     }
-    Opm::GravityColumnSolver<TransportModel> colsolver(model, *grid->c_grid());
+    Opm::GravityColumnSolver<TransportModel> colsolver(model, *grid->c_grid(), nltol, maxit);
 
     // State-related and source-related variables init.
     int num_cells = grid->c_grid()->number_of_cells;
@@ -418,20 +418,19 @@ main(int argc, char** argv)
 	    }
 	    break;
 	}
-    
     case 3:
 	{
-        std::cout << "==== Scenario 2: gravity convection.\n";
+	    std::cout << "==== Scenario 3: gravity segregation.\n";
 	    if (!use_gravity) {
-            std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
+		std::cout << "**** Warning: running gravity segregation scenario, but gravity is zero." << std::endl;
 	    }
 	    if (use_deck) {
-            std::cout << "**** Warning: running gravity convection scenario, which expects a cartesian grid."
+		std::cout << "**** Warning: running gravity segregation scenario, which expects a cartesian grid."
 			  << std::endl;
 	    }
 	    std::vector<double>& sat = state.saturation();
 	    const int *glob_cell = grid->c_grid()->global_cell;
-        // Heavy on top
+	    // Water on top
 	    for (int cell = 0; cell < num_cells; ++cell) {
 		const int* cd = grid->c_grid()->cartdims;
 		const int gc = glob_cell == 0 ? cell : glob_cell[cell];

opm/core/transport/GravityColumnSolver.hpp

@@ -36,7 +36,9 @@ namespace Opm
 	/// Note: the model will be changed since it stores computed
 	/// quantities in itself, such as mobilities.
 	GravityColumnSolver(Model& model,
-			    const UnstructuredGrid& grid);
+			    const UnstructuredGrid& grid,
+			    const double tol,
+			    const int maxit);
 
 	/// \param[in] columns         for each column (with logical cartesian indices as key),
 	///                            contains the cells on which to solve the segregation
@@ -51,10 +53,11 @@ namespace Opm
 			       const double dt,
 			       std::vector<double>& s,
 			       std::vector<double>& sol_vec);
-
 	Model& model_;
 	const UnstructuredGrid& grid_;
-    };
+	const double tol_;
+	const int maxit_;
+};
 
 } // namespace Opm
 

opm/core/transport/GravityColumnSolver_impl.hpp

@@ -26,8 +26,10 @@ namespace Opm
 
     template <class Model>
     GravityColumnSolver<Model>::GravityColumnSolver(Model& model,
-						    const UnstructuredGrid& grid)
+						    const UnstructuredGrid& grid,
-	: model_(model), grid_(grid)
+						    const double tol,
+						    const int maxit)
+	: model_(model), grid_(grid), tol_(tol), maxit_(maxit)
     {
     }
 
@@ -77,28 +79,30 @@ namespace Opm
 	// Initialize model. These things are done for the whole grid!
 	StateWithZeroFlux state(s); // This holds s by reference.
 	JacSys sys(grid_.number_of_cells);
+	std::vector<double> increment(grid_.number_of_cells, 0.0);
 	model_.initStep(state, grid_, sys);
 
-	const int max_iter = 40;
-	const double tol = 1e-4;
 	int iter = 0;
 	double max_delta = 1e100;
-	while (iter < max_iter) {
+	while (iter < maxit_) {
 	    model_.initIteration(state, grid_, sys);
 	    std::map<int, std::vector<int> >::const_iterator it;
 	    for (it = columns.begin(); it != columns.end(); ++it) {
-		solveSingleColumn(it->second, dt, s, sys.vector().writableSolution());
+		solveSingleColumn(it->second, dt, s, increment);
 	    }
-	    const double maxelem = *std::max_element(sys.vector().solution().begin(), sys.vector().solution().end());
+	    for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
-	    const double minelem = *std::min_element(sys.vector().solution().begin(), sys.vector().solution().end());
+		sys.vector().writableSolution()[cell] += increment[cell];
+	    }
+	    const double maxelem = *std::max_element(increment.begin(), increment.end());
+	    const double minelem = *std::min_element(increment.begin(), increment.end());
 	    max_delta = std::max(maxelem, -minelem);
 	    std::cout << "Iteration " << iter << "   max_delta = " << max_delta << std::endl;
-	    if (max_delta < tol) {
+	    if (max_delta < tol_) {
 		break;
 	    }
 	    ++iter;
 	}
-	if (max_delta >= tol) {
+	if (max_delta >= tol_) {
 	    THROW("Failed to converge!");
 	}
 	// Finalize.
@@ -168,7 +172,7 @@ namespace Opm
 	    THROW("Lapack reported error in dgtsv: " << info);
 	}
 	for (int ci = 0; ci < col_size; ++ci) {
-	    sol_vec[column_cells[ci]] = rhs[ci];
+	    sol_vec[column_cells[ci]] = -rhs[ci];
 	}
     }
 

opm/core/transport/reorder/TransportModelTwophase.cpp

@@ -28,7 +28,7 @@
 #include <numeric>
 
 
-// #define EXPERIMENT_GAUSS_SEIDEL
+#define EXPERIMENT_GAUSS_SEIDEL
 
 
 namespace Opm
@@ -302,17 +302,46 @@ namespace Opm
 	// Must store s0 before we start.
 	std::vector<double> s0(num_cells);
 	// Must set initial fractional flows before we start.
+	// Also, we compute the # of upstream neighbours.
+	// std::vector<int> num_upstream(num_cells);
 	for (int i = 0; i < num_cells; ++i) {
 	    const int cell = cells[i];
 	    fractionalflow_[cell] = fracFlow(saturation_[cell], cell);
 	    s0[i] = saturation_[cell];
+	    // num_upstream[i] = ia_upw_[cell + 1] - ia_upw_[cell];
 	}
 	// Solve once in each cell.
+	// std::vector<int> fully_marked_stack;
+	// fully_marked_stack.reserve(num_cells);
 	int num_iters = 0;
 	int update_count = 0; // Change name/meaning to cells_updated?
 	do {
 	    update_count = 0; // Must reset count for every iteration.
 	    for (int i = 0; i < num_cells; ++i) {
+		// while (!fully_marked_stack.empty()) {
+		//     // std::cout << "# fully marked cells = " << fully_marked_stack.size() << std::endl;
+		//     const int fully_marked_ci = fully_marked_stack.back();
+		//     fully_marked_stack.pop_back();
+		//     ++update_count;
+		//     const int cell = cells[fully_marked_ci];
+		//     const double old_s = saturation_[cell];
+		//     saturation_[cell] = s0[fully_marked_ci];
+		//     solveSingleCell(cell);
+		//     const double s_change = std::fabs(saturation_[cell] - old_s);
+		//     if (s_change > tol) {
+		// 	// Mark downwind cells.
+		// 	for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) {
+		// 	    const int downwind_cell = ja_downw_[j];
+		// 	    int ci = pos[downwind_cell];
+		// 	    ++needs_update[ci];
+		// 	    if (needs_update[ci] == num_upstream[ci]) {
+		// 		fully_marked_stack.push_back(ci);
+		// 	    }
+		// 	}
+		//     }
+		//     // Unmark this cell.
+		//     needs_update[fully_marked_ci] = 0;
+		// }
 		if (!needs_update[i]) {
 		    continue;
 		}
@@ -326,14 +355,23 @@ namespace Opm
 		    // Mark downwind cells.
 		    for (int j = ia_downw_[cell]; j < ia_downw_[cell+1]; ++j) {
 			const int downwind_cell = ja_downw_[j];
-			needs_update[pos[downwind_cell]] = 1;
+			int ci = pos[downwind_cell];
+			needs_update[ci] = 1;
+			// ++needs_update[ci];
+			// if (needs_update[ci] == num_upstream[ci]) {
+			//     fully_marked_stack.push_back(ci);
+			// }
 		    }
 		}
 		// Unmark this cell.
 		needs_update[i] = 0;
 	    }
-	    // std::cout << "Iter = " << num_iters << "    update_count = " << update_count << std::endl;
+	    // std::cout << "Iter = " << num_iters << "    update_count = " << update_count
+	    // 	      << "    # marked cells = "
+	    // 	      << std::accumulate(needs_update.begin(), needs_update.end(), 0) << std::endl;
 	} while (update_count > 0 && ++num_iters < max_iters);
+
+	// Done with iterations, check if we succeeded.
 	if (update_count > 0) {
 	    THROW("In solveMultiCell(), we did not converge after "
 	    	  << num_iters << " iterations. Remaining update count = " << update_count);