diff --git a/opm/core/transport/reorder/TransportModelInterface.cpp b/opm/core/transport/reorder/TransportModelInterface.cpp index cfaed289f..c52b528d6 100644 --- a/opm/core/transport/reorder/TransportModelInterface.cpp +++ b/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 sequence(grid.number_of_cells); std::vector 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. - assert(ncomponents == grid.number_of_cells); - for (int i = 0; i < grid.number_of_cells; ++i) { + // Invoke appropriate solve method for each interdependent component. + for (int comp = 0; comp < ncomponents; ++comp) { #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]]); + } } } diff --git a/opm/core/transport/reorder/TransportModelInterface.hpp b/opm/core/transport/reorder/TransportModelInterface.hpp index 6ad13aac6..4a7567095 100644 --- a/opm/core/transport/reorder/TransportModelInterface.hpp +++ b/opm/core/transport/reorder/TransportModelInterface.hpp @@ -26,17 +26,18 @@ namespace Opm { /// Interface for reordering transport models. - /// A transport model must provide the solveSingleCell() - /// method, and is expected to implement a solve() method - /// that will have an interface geared to the model's needs. - /// (The solve() method is therefore not virtual in this class). - /// The reorderAndTransport() method is provided as an - /// aid to implementing solve() in subclasses. + /// A transport model must provide the solveSingleCell() and + /// solveMultiCell methods, and is expected to implement a solve() + /// method that will have an interface geared to the model's + /// needs. (The solve() method is therefore not virtual in this + /// class.) The reorderAndTransport() method is provided as an aid + /// 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); }; diff --git a/opm/core/transport/reorder/TransportModelTwophase.cpp b/opm/core/transport/reorder/TransportModelTwophase.cpp index 2e56fae72..6a3216ff5 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.cpp +++ b/opm/core/transport/reorder/TransportModelTwophase.cpp @@ -23,6 +23,9 @@ #include #include +#include +#include + 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(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 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 }; diff --git a/opm/core/transport/reorder/TransportModelTwophase.hpp b/opm/core/transport/reorder/TransportModelTwophase.hpp index f1b88cd77..04b90bcdf 100644 --- a/opm/core/transport/reorder/TransportModelTwophase.hpp +++ b/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_;