OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Corrected various bugs in gravitation solver.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-03-20 17:43:25 +01:00
parent bd8288077f
commit 62a1c1ad98
2 changed files with 10 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
opm/polymer/GravityColumnSolverPolymer_impl.hpp
View File

@ -178,10 +178,11 @@ namespace Opm
const int kl = 3; const int kl = 3;
const int ku = 3; const int ku = 3;
const int nrow = 2*kl + ku + 1; const int nrow = 2*kl + ku + 1;
const BandMatrixCoeff bmc(col_size, ku, kl); const int N = 2*col_size; // N unknowns: s and c for each cell.
const BandMatrixCoeff bmc(N, ku, kl);
std::vector<double> hm(nrow*col_size, 0.0); // band matrix with 3 upper and 3 lower diagonals. std::vector<double> hm(nrow*N, 0.0); // band matrix with 3 upper and 3 lower diagonals.
std::vector<double> rhs(col_size, 0.0); std::vector<double> rhs(N, 0.0);
for (int ci = 0; ci < col_size; ++ci) { for (int ci = 0; ci < col_size; ++ci) {
std::vector<double> F(2, 0.); std::vector<double> F(2, 0.);
std::vector<double> dFd1(4, 0.); std::vector<double> dFd1(4, 0.);
@ -203,13 +204,13 @@ namespace Opm
if (c1 == prev_cell || c2 == prev_cell) { if (c1 == prev_cell || c2 == prev_cell) {
hm[bmc(2*ci + 0, 2*(ci - 1) + 0)] += dFd2[0]; hm[bmc(2*ci + 0, 2*(ci - 1) + 0)] += dFd2[0];
hm[bmc(2*ci + 0, 2*(ci - 1) + 1)] += dFd2[1]; hm[bmc(2*ci + 0, 2*(ci - 1) + 1)] += dFd2[1];
hm[bmc(2*ci + 1, 2*(ci - 1) + 1)] += dFd2[2]; hm[bmc(2*ci + 1, 2*(ci - 1) + 0)] += dFd2[2];
hm[bmc(2*ci + 1, 2*(ci - 1) + 1)] += dFd2[3]; hm[bmc(2*ci + 1, 2*(ci - 1) + 1)] += dFd2[3];
} else { } else {
ASSERT(c1 == next_cell || c2 == next_cell); ASSERT(c1 == next_cell || c2 == next_cell);
hm[bmc(2*ci + 0, 2*(ci + 1) + 0)] += dFd2[0]; hm[bmc(2*ci + 0, 2*(ci + 1) + 0)] += dFd2[0];
hm[bmc(2*ci + 0, 2*(ci + 1) + 1)] += dFd2[1]; hm[bmc(2*ci + 0, 2*(ci + 1) + 1)] += dFd2[1];
hm[bmc(2*ci + 1, 2*(ci + 1) + 1)] += dFd2[2]; hm[bmc(2*ci + 1, 2*(ci + 1) + 0)] += dFd2[2];
hm[bmc(2*ci + 1, 2*(ci + 1) + 1)] += dFd2[3]; hm[bmc(2*ci + 1, 2*(ci + 1) + 1)] += dFd2[3];
} }
hm[bmc(2*ci + 0, 2*ci + 0)] += dFd1[0]; hm[bmc(2*ci + 0, 2*ci + 0)] += dFd1[0];
@ -229,7 +230,7 @@ namespace Opm
hm[bmc(2*ci + 1, 2*ci + 0)] += dF[2]; hm[bmc(2*ci + 1, 2*ci + 0)] += dF[2];
hm[bmc(2*ci + 1, 2*ci + 1)] += dF[3]; hm[bmc(2*ci + 1, 2*ci + 1)] += dF[3];
rhs[2*ci + 0] += F[0]; rhs[2*ci + 0] += F[0];
rhs[2*ci + 1] += F[1]; rhs[2*ci + 1] += F[1];
} }
@ -239,7 +240,7 @@ namespace Opm
int info = 0; int info = 0;
int ipiv; int ipiv;
// Solution will be written to rhs. // Solution will be written to rhs.
dgbsv_(&col_size, &kl, &ku, &num_rhs, &hm[0], &nrow, &ipiv, &rhs[0], &col_size, &info); dgbsv_(&N, &kl, &ku, &num_rhs, &hm[0], &nrow, &ipiv, &rhs[0], &N, &info);
if (info != 0) { if (info != 0) {
THROW("Lapack reported error in dgtsv: " << info); THROW("Lapack reported error in dgtsv: " << info);
} }

5
opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp
View File

@ -441,7 +441,7 @@ namespace Opm {
typename JacobianSystem::vector_type& x = typename JacobianSystem::vector_type& x =
sys.vector().writableSolution(); sys.vector().writableSolution();
assert (x.size() == (::std::size_t) (g.number_of_cells)); assert (x.size() == (::std::size_t) (2*g.number_of_cells));
if (init_step_use_previous_sol_) { if (init_step_use_previous_sol_) {
std::fill(x.begin(), x.end(), 0.0); std::fill(x.begin(), x.end(), 0.0);
@ -546,10 +546,9 @@ namespace Opm {
double *s = &state.saturation()[0*2 + 0]; double *s = &state.saturation()[0*2 + 0];
double *c = &state.concentration()[0*1 + 0]; double *c = &state.concentration()[0*1 + 0];
for (int cell = 0; cell < g.number_of_cells; ++cell, s += 2) { for (int cell = 0; cell < g.number_of_cells; ++cell, s += 2, c += 1) {
s[0] += x[2*cell + 0]; s[0] += x[2*cell + 0];
c[0] += x[2*cell + 1]; c[0] += x[2*cell + 1];
c += 1;
double s_min = fluid_.s_min(cell); double s_min = fluid_.s_min(cell);
double s_max = fluid_.s_max(cell); double s_max = fluid_.s_max(cell);
assert(s[0] >= s_min - sat_tol_); assert(s[0] >= s_min - sat_tol_);