00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include <opm/polymer/GravityColumnSolverPolymer.hpp>
00040 #include <opm/core/linalg/blas_lapack.h>
00041 #include <opm/common/ErrorMacros.hpp>
00042 #include <iterator>
00043 #include <iostream>
00044 #include <cmath>
00045 #include <algorithm>
00046
00047 namespace Opm
00048 {
00049
00050 template <class FluxModel, class Model>
00051 GravityColumnSolverPolymer<FluxModel, Model>::GravityColumnSolverPolymer(FluxModel& fmodel,
00052 const Model& model,
00053 const UnstructuredGrid& grid,
00054 const double tol,
00055 const int maxit)
00056 : fmodel_(fmodel), model_(model), grid_(grid), tol_(tol), maxit_(maxit)
00057 {
00058 }
00059
00060 namespace {
00061 struct ZeroVec
00062 {
00063 double operator[](int) const { return 0.0; }
00064 };
00065 struct StateWithZeroFlux
00066 {
00067 StateWithZeroFlux(std::vector<double>& s, std::vector<double>& c, std::vector<double>& cmax_arg) : sat(s), cpoly(c), cmax(cmax_arg) {}
00068 const ZeroVec& faceflux() const { return zv; }
00069 const std::vector<double>& saturation() const { return sat; }
00070 std::vector<double>& saturation() { return sat; }
00071 const std::vector<double>& concentration() const { return cpoly; }
00072 std::vector<double>& concentration() { return cpoly; }
00073 const std::vector<double>& maxconcentration() const { return cmax; }
00074 std::vector<double>& maxconcentration() { return cmax; }
00075 ZeroVec zv;
00076 std::vector<double>& sat;
00077 std::vector<double>& cpoly;
00078 std::vector<double>& cmax;
00079 };
00080
00081 struct Vecs
00082 {
00083 Vecs(int sz) : sol(sz, 0.0) {}
00084 const std::vector<double>& solution() const { return sol; }
00085 std::vector<double>& writableSolution() { return sol; }
00086 std::vector<double> sol;
00087 };
00088 struct JacSys
00089 {
00090 JacSys(int sz) : v(sz) {}
00091 const Vecs& vector() const { return v; }
00092 Vecs& vector() { return v; }
00093 Vecs v;
00094 typedef std::vector<double> vector_type;
00095 };
00096
00097 struct BandMatrixCoeff
00098 {
00099 BandMatrixCoeff(int N, int ku, int kl) : N_(N), ku_(ku), kl_(kl), nrow_(2*kl + ku + 1) {
00100 }
00101
00102
00103
00104
00105 int operator ()(int i, int j) const {
00106 return kl_ + ku_ + i - j + j*nrow_;
00107 }
00108
00109 const int N_;
00110 const int ku_;
00111 const int kl_;
00112 const int nrow_;
00113 };
00114
00115 }
00116
00117
00118
00124 template <class FluxModel, class Model>
00125 void GravityColumnSolverPolymer<FluxModel, Model>::solve(const std::vector<std::vector<int> >& columns,
00126 const double dt,
00127 std::vector<double>& s,
00128 std::vector<double>& c,
00129 std::vector<double>& cmax
00130 )
00131 {
00132
00133 StateWithZeroFlux state(s, c, cmax);
00134 JacSys sys(2*grid_.number_of_cells);
00135 std::vector<double> increment(2*grid_.number_of_cells, 0.0);
00136 fmodel_.initStep(state, grid_, sys);
00137
00138 int iter = 0;
00139 double max_delta = 1e100;
00140 const double cmax_cell = 2.0*model_.cMax();
00141 const double tol_c_cell = 1e-2*cmax_cell;
00142 while (iter < maxit_) {
00143 fmodel_.initIteration(state, grid_, sys);
00144 int size = columns.size();
00145 for(int i = 0; i < size; ++i) {
00146 solveSingleColumn(columns[i], dt, s, c, cmax, increment);
00147 }
00148 for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
00149 double& s_cell = sys.vector().writableSolution()[2*cell + 0];
00150 double& c_cell = sys.vector().writableSolution()[2*cell + 1];
00151 s_cell += increment[2*cell + 0];
00152 c_cell += increment[2*cell + 1];
00153 if (s_cell < 0.) {
00154 double& incr = increment[2*cell + 0];
00155 s_cell -= incr;
00156 if (std::fabs(incr) < 1e-2) {
00157 incr = -s_cell;
00158 s_cell = 0.;
00159 } else {
00160 incr = -s_cell/2.0;
00161 s_cell = s_cell/2.0;
00162 }
00163 }
00164 if (s_cell > 1.) {
00165 double& incr = increment[2*cell + 0];
00166 s_cell -= incr;
00167 if (std::fabs(incr) < 1e-2) {
00168 incr = 1. - s_cell;
00169 s_cell = 1.;
00170 } else {
00171 incr = (1 - s_cell)/2.0;
00172 s_cell = (1 + s_cell)/2.0;
00173 }
00174 }
00175 if (c_cell < 0.) {
00176 double& incr = increment[2*cell + 1];
00177 c_cell -= incr;
00178 if (std::fabs(incr) < tol_c_cell) {
00179 incr = -c_cell;
00180 c_cell = 0.;
00181 } else {
00182 incr = -c_cell/2.0;
00183 c_cell = c_cell/2.0;
00184 }
00185 }
00186 if (c_cell > cmax_cell) {
00187 double& incr = increment[2*cell + 1];
00188 c_cell -= incr;
00189 if (std::fabs(incr) < tol_c_cell) {
00190 incr = cmax_cell - c_cell;
00191 c_cell = cmax_cell;
00192 } else {
00193 incr = (cmax_cell - c_cell)/2.0;
00194 c_cell = (cmax_cell + c_cell)/2.0;
00195 }
00196 }
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212 }
00213 const double maxelem = *std::max_element(increment.begin(), increment.end());
00214 const double minelem = *std::min_element(increment.begin(), increment.end());
00215 max_delta = std::max(maxelem, -minelem);
00216 std::cout << "Iteration " << iter << " max_delta = " << max_delta << std::endl;
00217 if (max_delta < tol_) {
00218 break;
00219 }
00220 ++iter;
00221 }
00222 if (max_delta >= tol_) {
00223 OPM_THROW(std::runtime_error, "Failed to converge!");
00224 }
00225
00226
00227
00228
00229 fmodel_.finishStep(grid_, sys.vector().solution(), state);
00230 }
00231
00232
00233
00234
00238 template <class FluxModel, class Model>
00239 void GravityColumnSolverPolymer<FluxModel, Model>::solveSingleColumn(const std::vector<int>& column_cells,
00240 const double dt,
00241 std::vector<double>& s,
00242 std::vector<double>& c,
00243 std::vector<double>& cmax,
00244 std::vector<double>& sol_vec)
00245 {
00246
00247
00248 int col_size = column_cells.size();
00249
00250
00251
00252
00253
00254
00255
00256 StateWithZeroFlux state(s, c, cmax);
00257
00258
00259 const int kl = 3;
00260 const int ku = 3;
00261 const int nrow = 2*kl + ku + 1;
00262 const int N = 2*col_size;
00263 std::vector<double> hm(nrow*N, 0.0);
00264 std::vector<double> rhs(N, 0.0);
00265 const BandMatrixCoeff bmc(N, ku, kl);
00266
00267
00268 for (int ci = 0; ci < col_size; ++ci) {
00269 std::vector<double> F(2, 0.);
00270 std::vector<double> dFd1(4, 0.);
00271 std::vector<double> dFd2(4, 0.);
00272 std::vector<double> dF(4, 0.);
00273 const int cell = column_cells[ci];
00274 const int prev_cell = (ci == 0) ? -999 : column_cells[ci - 1];
00275 const int next_cell = (ci == col_size - 1) ? -999 : column_cells[ci + 1];
00276
00277 for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) {
00278 const int face = grid_.cell_faces[j];
00279 const int c1 = grid_.face_cells[2*face + 0];
00280 const int c2 = grid_.face_cells[2*face + 1];
00281 if (c1 == prev_cell || c2 == prev_cell || c1 == next_cell || c2 == next_cell) {
00282 F.assign(2, 0.);
00283 dFd1.assign(4, 0.);
00284 dFd2.assign(4, 0.);
00285 fmodel_.fluxConnection(state, grid_, dt, cell, face, &F[0], &dFd1[0], &dFd2[0]);
00286 if (c1 == prev_cell || c2 == prev_cell) {
00287 hm[bmc(2*ci + 0, 2*(ci - 1) + 0)] += dFd2[0];
00288 hm[bmc(2*ci + 0, 2*(ci - 1) + 1)] += dFd2[1];
00289 hm[bmc(2*ci + 1, 2*(ci - 1) + 0)] += dFd2[2];
00290 hm[bmc(2*ci + 1, 2*(ci - 1) + 1)] += dFd2[3];
00291 } else {
00292 assert(c1 == next_cell || c2 == next_cell);
00293 hm[bmc(2*ci + 0, 2*(ci + 1) + 0)] += dFd2[0];
00294 hm[bmc(2*ci + 0, 2*(ci + 1) + 1)] += dFd2[1];
00295 hm[bmc(2*ci + 1, 2*(ci + 1) + 0)] += dFd2[2];
00296 hm[bmc(2*ci + 1, 2*(ci + 1) + 1)] += dFd2[3];
00297 }
00298 hm[bmc(2*ci + 0, 2*ci + 0)] += dFd1[0];
00299 hm[bmc(2*ci + 0, 2*ci + 1)] += dFd1[1];
00300 hm[bmc(2*ci + 1, 2*ci + 0)] += dFd1[2];
00301 hm[bmc(2*ci + 1, 2*ci + 1)] += dFd1[3];
00302
00303 rhs[2*ci + 0] += F[0];
00304 rhs[2*ci + 1] += F[1];
00305 }
00306 }
00307 F.assign(2, 0.);
00308 dF.assign(4, 0.);
00309 fmodel_.accumulation(grid_, cell, &F[0], &dF[0]);
00310 hm[bmc(2*ci + 0, 2*ci + 0)] += dF[0];
00311 hm[bmc(2*ci + 0, 2*ci + 1)] += dF[1];
00312 hm[bmc(2*ci + 1, 2*ci + 0)] += dF[2];
00313 if (std::abs(dF[3]) < 1e-12) {
00314 hm[bmc(2*ci + 1, 2*ci + 1)] += 1e-12;
00315 } else {
00316 hm[bmc(2*ci + 1, 2*ci + 1)] += dF[3];
00317 }
00318
00319 rhs[2*ci + 0] += F[0];
00320 rhs[2*ci + 1] += F[1];
00321
00322 }
00323
00324
00325 const int num_rhs = 1;
00326 int info = 0;
00327 std::vector<int> ipiv(N, 0);
00328
00329 dgbsv_(&N, &kl, &ku, &num_rhs, &hm[0], &nrow, &ipiv[0], &rhs[0], &N, &info);
00330 if (info != 0) {
00331 std::cerr << "Failed column cells: ";
00332 std::copy(column_cells.begin(), column_cells.end(), std::ostream_iterator<int>(std::cerr, " "));
00333 std::cerr << "\n";
00334 OPM_THROW(std::runtime_error, "Lapack reported error in dgtsv: " << info);
00335 }
00336 for (int ci = 0; ci < col_size; ++ci) {
00337 sol_vec[2*column_cells[ci] + 0] = -rhs[2*ci + 0];
00338 sol_vec[2*column_cells[ci] + 1] = -rhs[2*ci + 1];
00339 }
00340 }
00341
00342 }