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
47namespace 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
N_
const int N_
Definition: GravityColumnSolverPolymer_impl.hpp:109
nrow_
const int nrow_
Definition: GravityColumnSolverPolymer_impl.hpp:112
sol
std::vector< double > sol
Definition: GravityColumnSolverPolymer_impl.hpp:86
zv
ZeroVec zv
Definition: GravityColumnSolverPolymer_impl.hpp:75
kl_
const int kl_
Definition: GravityColumnSolverPolymer_impl.hpp:111
cpoly
std::vector< double > & cpoly
Definition: GravityColumnSolverPolymer_impl.hpp:77
sat
std::vector< double > & sat
Definition: GravityColumnSolverPolymer_impl.hpp:76
ku_
const int ku_
Definition: GravityColumnSolverPolymer_impl.hpp:110
v
Vecs v
Definition: GravityColumnSolverPolymer_impl.hpp:93
cmax
std::vector< double > & cmax
Definition: GravityColumnSolverPolymer_impl.hpp:78
Opm::GravityColumnSolverPolymer
Definition: GravityColumnSolverPolymer.hpp:34
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
Opm::GravityColumnSolverPolymer::GravityColumnSolverPolymer
GravityColumnSolverPolymer(FluxModel &fmodel, const Model &model, const UnstructuredGrid &grid, const double tol, const int maxit)
Definition: GravityColumnSolverPolymer_impl.hpp:51
Opm
Definition: CompressibleTpfaPolymer.hpp:33