GravityColumnSolverPolymer_impl.hpp
Go to the documentation of this file.
1 /*
2  Copyright 2012 SINTEF ICT, Applied Mathematics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 /*
21  Copyright 2012 SINTEF ICT, Applied Mathematics.
22 
23  This file is part of the Open Porous Media project (OPM).
24 
25  OPM is free software: you can redistribute it and/or modify
26  it under the terms of the GNU General Public License as published by
27  the Free Software Foundation, either version 3 of the License, or
28  (at your option) any later version.
29 
30  OPM is distributed in the hope that it will be useful,
31  but WITHOUT ANY WARRANTY; without even the implied warranty of
32  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33  GNU General Public License for more details.
34 
35  You should have received a copy of the GNU General Public License
36  along with OPM. If not, see <http://www.gnu.org/licenses/>.
37 */
38 
39 #include <opm/polymer/GravityColumnSolverPolymer.hpp>
40 #include <opm/core/linalg/blas_lapack.h>
41 #include <opm/common/ErrorMacros.hpp>
42 #include <iterator>
43 #include <iostream>
44 #include <cmath>
45 #include <algorithm>
46 
47 namespace Opm
48 {
49 
50  template <class FluxModel, class Model>
51  GravityColumnSolverPolymer<FluxModel, Model>::GravityColumnSolverPolymer(FluxModel& fmodel,
52  const Model& model,
53  const UnstructuredGrid& grid,
54  const double tol,
55  const int maxit)
56  : fmodel_(fmodel), model_(model), grid_(grid), tol_(tol), maxit_(maxit)
57  {
58  }
59 
60  namespace {
61  struct ZeroVec
62  {
63  double operator[](int) const { return 0.0; }
64  };
65  struct StateWithZeroFlux
66  {
67  StateWithZeroFlux(std::vector<double>& s, std::vector<double>& c, std::vector<double>& cmax_arg) : sat(s), cpoly(c), cmax(cmax_arg) {}
68  const ZeroVec& faceflux() const { return zv; }
69  const std::vector<double>& saturation() const { return sat; }
70  std::vector<double>& saturation() { return sat; }
71  const std::vector<double>& concentration() const { return cpoly; }
72  std::vector<double>& concentration() { return cpoly; }
73  const std::vector<double>& maxconcentration() const { return cmax; }
74  std::vector<double>& maxconcentration() { return cmax; }
75  ZeroVec zv;
76  std::vector<double>& sat;
77  std::vector<double>& cpoly;
78  std::vector<double>& cmax;
79  };
80 
81  struct Vecs
82  {
83  Vecs(int sz) : sol(sz, 0.0) {}
84  const std::vector<double>& solution() const { return sol; }
85  std::vector<double>& writableSolution() { return sol; }
86  std::vector<double> sol;
87  };
88  struct JacSys
89  {
90  JacSys(int sz) : v(sz) {}
91  const Vecs& vector() const { return v; }
92  Vecs& vector() { return v; }
93  Vecs v;
94  typedef std::vector<double> vector_type;
95  };
96 
97  struct BandMatrixCoeff
98  {
99  BandMatrixCoeff(int N, int ku, int kl) : N_(N), ku_(ku), kl_(kl), nrow_(2*kl + ku + 1) {
100  }
101 
102 
103  // compute the position where to store the coefficient of a matrix A_{i,j} (i,j=0,...,N-1)
104  // in a array which is sent to the band matrix solver of LAPACK.
105  int operator ()(int i, int j) const {
106  return kl_ + ku_ + i - j + j*nrow_;
107  }
108 
109  const int N_;
110  const int ku_;
111  const int kl_;
112  const int nrow_;
113  };
114 
115  } // anon namespace
116 
117 
118 
124  template <class FluxModel, class Model>
125  void GravityColumnSolverPolymer<FluxModel, Model>::solve(const std::vector<std::vector<int> >& columns,
126  const double dt,
127  std::vector<double>& s,
128  std::vector<double>& c,
129  std::vector<double>& cmax
130  )
131  {
132  // Initialize model. These things are done for the whole grid!
133  StateWithZeroFlux state(s, c, cmax); // This holds s, c and cmax by reference.
134  JacSys sys(2*grid_.number_of_cells);
135  std::vector<double> increment(2*grid_.number_of_cells, 0.0);
136  fmodel_.initStep(state, grid_, sys);
137 
138  int iter = 0;
139  double max_delta = 1e100;
140  const double cmax_cell = 2.0*model_.cMax();
141  const double tol_c_cell = 1e-2*cmax_cell;
142  while (iter < maxit_) {
143  fmodel_.initIteration(state, grid_, sys);
144  int size = columns.size();
145  for(int i = 0; i < size; ++i) {
146  solveSingleColumn(columns[i], dt, s, c, cmax, increment);
147  }
148  for (int cell = 0; cell < grid_.number_of_cells; ++cell) {
149  double& s_cell = sys.vector().writableSolution()[2*cell + 0];
150  double& c_cell = sys.vector().writableSolution()[2*cell + 1];
151  s_cell += increment[2*cell + 0];
152  c_cell += increment[2*cell + 1];
153  if (s_cell < 0.) {
154  double& incr = increment[2*cell + 0];
155  s_cell -= incr;
156  if (std::fabs(incr) < 1e-2) {
157  incr = -s_cell;
158  s_cell = 0.;
159  } else {
160  incr = -s_cell/2.0;
161  s_cell = s_cell/2.0;
162  }
163  }
164  if (s_cell > 1.) {
165  double& incr = increment[2*cell + 0];
166  s_cell -= incr;
167  if (std::fabs(incr) < 1e-2) {
168  incr = 1. - s_cell;
169  s_cell = 1.;
170  } else {
171  incr = (1 - s_cell)/2.0;
172  s_cell = (1 + s_cell)/2.0;
173  }
174  }
175  if (c_cell < 0.) {
176  double& incr = increment[2*cell + 1];
177  c_cell -= incr;
178  if (std::fabs(incr) < tol_c_cell) {
179  incr = -c_cell;
180  c_cell = 0.;
181  } else {
182  incr = -c_cell/2.0;
183  c_cell = c_cell/2.0;
184  }
185  }
186  if (c_cell > cmax_cell) {
187  double& incr = increment[2*cell + 1];
188  c_cell -= incr;
189  if (std::fabs(incr) < tol_c_cell) {
190  incr = cmax_cell - c_cell;
191  c_cell = cmax_cell;
192  } else {
193  incr = (cmax_cell - c_cell)/2.0;
194  c_cell = (cmax_cell + c_cell)/2.0;
195  }
196  }
197 
198  // if (s_cell < 0.) {
199  // increment[2*cell + 0] = increment[2*cell + 0] - s_cell;
200  // s_cell = 0.;
201  // } else if (s_cell > 1.) {
202  // increment[2*cell + 0] = increment[2*cell + 0] - s_cell + 1.;
203  // s_cell = 1.;
204  // }
205  // if (c_cell < 0) {
206  // increment[2*cell + 1] = increment[2*cell + 1] - c_cell;
207  // c_cell = 0.;
208  // } else if (c_cell > cmax_cell) {
209  // increment[2*cell + 1] = increment[2*cell + 1] - c_cell + cmax_cell;
210  // c_cell = cmax_cell;
211  // }
212  }
213  const double maxelem = *std::max_element(increment.begin(), increment.end());
214  const double minelem = *std::min_element(increment.begin(), increment.end());
215  max_delta = std::max(maxelem, -minelem);
216  std::cout << "Iteration " << iter << " max_delta = " << max_delta << std::endl;
217  if (max_delta < tol_) {
218  break;
219  }
220  ++iter;
221  }
222  if (max_delta >= tol_) {
223  OPM_THROW(std::runtime_error, "Failed to converge!");
224  }
225  // Finalize.
226  // fmodel_.finishIteration(); //
227  // finishStep() writes to state, which holds s by reference.
228  // This will update the entire grid's state... cmax is updated here.
229  fmodel_.finishStep(grid_, sys.vector().solution(), state);
230  }
231 
232 
233 
234 
238  template <class FluxModel, class Model>
239  void GravityColumnSolverPolymer<FluxModel, Model>::solveSingleColumn(const std::vector<int>& column_cells,
240  const double dt,
241  std::vector<double>& s,
242  std::vector<double>& c,
243  std::vector<double>& cmax,
244  std::vector<double>& sol_vec)
245  {
246  // This is written only to work with SinglePointUpwindTwoPhase,
247  // not with arbitrary problem models.
248  int col_size = column_cells.size();
249 
250  // if (col_size == 1) {
251  // sol_vec[2*column_cells[0] + 0] = 0.0;
252  // sol_vec[2*column_cells[0] + 1] = 0.0;
253  // return;
254  // }
255 
256  StateWithZeroFlux state(s, c, cmax); // This holds s by reference.
257 
258  // Assemble.
259  const int kl = 3;
260  const int ku = 3;
261  const int nrow = 2*kl + ku + 1;
262  const int N = 2*col_size; // N unknowns: s and c for each cell.
263  std::vector<double> hm(nrow*N, 0.0); // band matrix with 3 upper and 3 lower diagonals.
264  std::vector<double> rhs(N, 0.0);
265  const BandMatrixCoeff bmc(N, ku, kl);
266 
267 
268  for (int ci = 0; ci < col_size; ++ci) {
269  std::vector<double> F(2, 0.);
270  std::vector<double> dFd1(4, 0.);
271  std::vector<double> dFd2(4, 0.);
272  std::vector<double> dF(4, 0.);
273  const int cell = column_cells[ci];
274  const int prev_cell = (ci == 0) ? -999 : column_cells[ci - 1];
275  const int next_cell = (ci == col_size - 1) ? -999 : column_cells[ci + 1];
276  // model_.initResidual(cell, F);
277  for (int j = grid_.cell_facepos[cell]; j < grid_.cell_facepos[cell+1]; ++j) {
278  const int face = grid_.cell_faces[j];
279  const int c1 = grid_.face_cells[2*face + 0];
280  const int c2 = grid_.face_cells[2*face + 1];
281  if (c1 == prev_cell || c2 == prev_cell || c1 == next_cell || c2 == next_cell) {
282  F.assign(2, 0.);
283  dFd1.assign(4, 0.);
284  dFd2.assign(4, 0.);
285  fmodel_.fluxConnection(state, grid_, dt, cell, face, &F[0], &dFd1[0], &dFd2[0]);
286  if (c1 == prev_cell || c2 == prev_cell) {
287  hm[bmc(2*ci + 0, 2*(ci - 1) + 0)] += dFd2[0];
288  hm[bmc(2*ci + 0, 2*(ci - 1) + 1)] += dFd2[1];
289  hm[bmc(2*ci + 1, 2*(ci - 1) + 0)] += dFd2[2];
290  hm[bmc(2*ci + 1, 2*(ci - 1) + 1)] += dFd2[3];
291  } else {
292  assert(c1 == next_cell || c2 == next_cell);
293  hm[bmc(2*ci + 0, 2*(ci + 1) + 0)] += dFd2[0];
294  hm[bmc(2*ci + 0, 2*(ci + 1) + 1)] += dFd2[1];
295  hm[bmc(2*ci + 1, 2*(ci + 1) + 0)] += dFd2[2];
296  hm[bmc(2*ci + 1, 2*(ci + 1) + 1)] += dFd2[3];
297  }
298  hm[bmc(2*ci + 0, 2*ci + 0)] += dFd1[0];
299  hm[bmc(2*ci + 0, 2*ci + 1)] += dFd1[1];
300  hm[bmc(2*ci + 1, 2*ci + 0)] += dFd1[2];
301  hm[bmc(2*ci + 1, 2*ci + 1)] += dFd1[3];
302 
303  rhs[2*ci + 0] += F[0];
304  rhs[2*ci + 1] += F[1];
305  }
306  }
307  F.assign(2, 0.);
308  dF.assign(4, 0.);
309  fmodel_.accumulation(grid_, cell, &F[0], &dF[0]);
310  hm[bmc(2*ci + 0, 2*ci + 0)] += dF[0];
311  hm[bmc(2*ci + 0, 2*ci + 1)] += dF[1];
312  hm[bmc(2*ci + 1, 2*ci + 0)] += dF[2];
313  if (std::abs(dF[3]) < 1e-12) {
314  hm[bmc(2*ci + 1, 2*ci + 1)] += 1e-12;
315  } else {
316  hm[bmc(2*ci + 1, 2*ci + 1)] += dF[3];
317  }
318 
319  rhs[2*ci + 0] += F[0];
320  rhs[2*ci + 1] += F[1];
321 
322  }
323  // model_.sourceTerms(); // Not needed
324  // Solve.
325  const int num_rhs = 1;
326  int info = 0;
327  std::vector<int> ipiv(N, 0);
328  // Solution will be written to rhs.
329  dgbsv_(&N, &kl, &ku, &num_rhs, &hm[0], &nrow, &ipiv[0], &rhs[0], &N, &info);
330  if (info != 0) {
331  std::cerr << "Failed column cells: ";
332  std::copy(column_cells.begin(), column_cells.end(), std::ostream_iterator<int>(std::cerr, " "));
333  std::cerr << "\n";
334  OPM_THROW(std::runtime_error, "Lapack reported error in dgtsv: " << info);
335  }
336  for (int ci = 0; ci < col_size; ++ci) {
337  sol_vec[2*column_cells[ci] + 0] = -rhs[2*ci + 0];
338  sol_vec[2*column_cells[ci] + 1] = -rhs[2*ci + 1];
339  }
340  }
341 
342 } // namespace Opm
Opm
Definition: CompressibleTpfaPolymer.hpp:32
sat
std::vector< double > & sat
Definition: GravityColumnSolverPolymer_impl.hpp:76
N_
const int N_
Definition: GravityColumnSolverPolymer_impl.hpp:109
cmax
std::vector< double > & cmax
Definition: GravityColumnSolverPolymer_impl.hpp:78
GravityColumnSolverPolymer.hpp
kl_
const int kl_
Definition: GravityColumnSolverPolymer_impl.hpp:111
Opm::GravityColumnSolverPolymer::GravityColumnSolverPolymer
GravityColumnSolverPolymer(FluxModel &fmodel, const Model &model, const UnstructuredGrid &grid, const double tol, const int maxit)
Definition: GravityColumnSolverPolymer_impl.hpp:51
nrow_
const int nrow_
Definition: GravityColumnSolverPolymer_impl.hpp:112
cpoly
std::vector< double > & cpoly
Definition: GravityColumnSolverPolymer_impl.hpp:77
v
Vecs v
Definition: GravityColumnSolverPolymer_impl.hpp:93
ku_
const int ku_
Definition: GravityColumnSolverPolymer_impl.hpp:110
Opm::GravityColumnSolverPolymer
Definition: GravityColumnSolverPolymer.hpp:33
zv
ZeroVec zv
Definition: GravityColumnSolverPolymer_impl.hpp:75
Opm::GravityColumnSolverPolymer::solve
void solve(const std::vector< std::vector< int > > &columns, const double dt, std::vector< double > &s, std::vector< double > &c, std::vector< double > &cmax)
Definition: GravityColumnSolverPolymer_impl.hpp:125
sol
std::vector< double > sol
Definition: GravityColumnSolverPolymer_impl.hpp:86