From ccdacb791c84445ee47d1dbee4c43fed285f06ac Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Fri, 17 Feb 2012 09:39:45 +0100
Subject: [PATCH] Added lots of (inactive) experimental code. The code attempts
 to improve #iterations in the Gauss-Seidel-like multicell solver by improving
 ordering. In general, experiment failed to improve #iterations, except for
 one: totally random order was the best (for the 100x100 case tried)!

---
 .../reorder/TransportModelTwophase.cpp        | 80 ++++++++++++++++++-
 .../reorder/TransportModelTwophase.hpp        |  4 +
 2 files changed, 82 insertions(+), 2 deletions(-)

diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp
index 6a3216ff..80f704a7 100644
--- a/opm/core/transport/reorder/TransportModelTwophase.cpp
+++ b/opm/core/transport/reorder/TransportModelTwophase.cpp
@@ -20,7 +20,7 @@
 #include <opm/core/transport/reorder/TransportModelTwophase.hpp>
 #include <opm/core/fluid/IncompPropertiesInterface.hpp>
 #include <opm/core/grid.h>
-#include <opm/core/transport/reorder/nlsolvers.h>
+#include <opm/core/transport/reorder/reordersequence.h>
 #include <opm/core/utility/RootFinders.hpp>
 
 #include <fstream>
@@ -41,6 +41,8 @@ namespace Opm
 	  dt_(0.0),
 	  saturation_(0),
 	  fractionalflow_(grid.number_of_cells, -1.0)
+	  // ia_(grid.number_of_cells + 1, -1),
+	  // ja_(grid.number_of_faces, -1)
     {
 	if (props.numPhases() != 2) {
 	    THROW("Property object must have 2 phases");
@@ -57,6 +59,15 @@ namespace Opm
 	source_ = source;
 	dt_ = dt;
 	saturation_ = saturation;
+
+	// std::vector<int> seq(grid_.number_of_cells);
+	// std::vector<int> comp(grid_.number_of_cells + 1);
+	// int ncomp;
+	// compute_sequence_graph(&grid_, darcyflux_,
+	// 		       &seq[0], &comp[0], &ncomp,
+	// 		       &ia_[0], &ja_[0]);
+
+
 	reorderAndTransport(grid_, darcyflux);
     }
 
@@ -135,6 +146,71 @@ namespace Opm
     {
 	// std::ofstream os("dump");
 	// std::copy(cells, cells + num_cells, std::ostream_iterator<double>(os, "\n"));
+
+	// Experiment: try a breath-first search to build a more suitable ordering.
+	// Verdict: failed to improve #iterations.
+	// {
+	//     std::vector<int> pos(grid_.number_of_cells, -1);
+	//     for (int i = 0; i < num_cells; ++i) {
+	// 	const int cell = cells[i];
+	// 	pos[cell] = i;
+	//     }
+	//     std::vector<int> done_pos(num_cells, 0);
+	//     std::vector<int> upstream_pos;
+	//     std::vector<int> new_pos;
+	//     upstream_pos.push_back(0);
+	//     done_pos[0] = 1;
+	//     int current = 0;
+	//     while (new_pos.size() < num_cells) {
+	// 	const int i = upstream_pos[current++];
+	// 	new_pos.push_back(i);
+	// 	const int cell = cells[i];
+	// 	for (int j = ia_[cell]; j < ia_[cell+1]; ++j) {
+	// 	    const int opos = pos[ja_[j]];
+	// 	    if (!done_pos[opos]) {
+	// 		upstream_pos.push_back(opos);
+	// 		done_pos[opos] = 1;
+	// 	    }
+	// 	}
+	//     }
+	//     std::reverse(new_pos.begin(), new_pos.end());
+	//     std::copy(new_pos.begin(), new_pos.end(), const_cast<int*>(cells));
+	// }
+
+	// Experiment: try a random ordering.
+	// Verdict: amazingly, reduced #iterations by approx. 25%!
+	// int* c = const_cast<int*>(cells);
+	// std::random_shuffle(c, c + num_cells);
+
+	// Experiment: implement a metric measuring badness of ordering
+	//             as average distance in (cyclic) ordering from
+	//             upstream neighbours.
+	// Verdict: does not seem to predict #iterations very well, if at all.
+	// std::vector<int> pos(grid_.number_of_cells, -1);
+	// for (int i = 0; i < num_cells; ++i) {
+	//     const int cell = cells[i];
+	//     pos[cell] = i;
+	// }
+	// double diffsum = 0;
+	// for (int i = 0; i < num_cells; ++i) {
+	//     const int cell = cells[i];
+	//     int num_upstream = 0;
+	//     int loc_diffsum = 0;
+	//     for (int j = ia_[cell]; j < ia_[cell+1]; ++j) {
+	// 	const int opos = pos[ja_[j]];
+	// 	if (opos == -1) {
+	// 	    std::cout << "Hmmm" << std::endl;
+	// 	    continue;
+	// 	}
+	// 	++num_upstream;
+	// 	const int diff = (i - opos + num_cells) % num_cells;
+	// 	loc_diffsum += diff;
+	//     }
+	//     diffsum += double(loc_diffsum)/double(num_upstream);
+	// }
+	// std::cout << "Average distance from upstream neighbours: " << diffsum/double(num_cells)
+	// 	  << std::endl;
+
 	double max_s_change = 0.0;
 	const double tol = 1e-9;
 	const int max_iters = 300;
@@ -155,7 +231,7 @@ namespace Opm
 		const double old_s = saturation_[cell];
 		saturation_[cell] = s0[i];
 		solveSingleCell(cell);
-		// std::cout << "delta s = " << saturation_[cell] - old_s << std::endl;
+		// std::cout << "cell = " << cell << "    delta s = " << saturation_[cell] - old_s << std::endl;
 		if (max_s_change < std::fabs(saturation_[cell] - old_s)) {
 		    max_change_cell = cell;
 		}
diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp
index 04b90bcd..c5ec589e 100644
--- a/opm/core/transport/reorder/TransportModelTwophase.hpp
+++ b/opm/core/transport/reorder/TransportModelTwophase.hpp
@@ -57,6 +57,10 @@ namespace Opm
 	std::vector<double> fractionalflow_;  // one per cell
 	const double* visc_;
 
+	// Storing the upwind graph for experiments.
+	// std::vector<int> ia_;
+	// std::vector<int> ja_;
+
 	struct Residual;
 	double fracFlow(double s, int cell) const;
     };