OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #73 from atgeirr/velocity_interpolation

Browse Source 
Fix output in case of LAPACK error.
This commit is contained in:
Bård Skaflestad 2012-10-17 03:53:27 -07:00
commit 73a836de8f
3 changed files with 7 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
View File

@ -465,6 +465,7 @@ namespace Opm
tof_coeff_ = &tof_coeff[0]; tof_coeff_ = &tof_coeff[0];
rhs_.resize(num_basis); rhs_.resize(num_basis);
jac_.resize(num_basis*num_basis); jac_.resize(num_basis*num_basis);
orig_jac_.resize(num_basis*num_basis);
basis_.resize(num_basis); basis_.resize(num_basis);
basis_nb_.resize(num_basis); basis_nb_.resize(num_basis);
grad_basis_.resize(num_basis*grid_.dimensions); grad_basis_.resize(num_basis*grid_.dimensions);
@ -617,6 +618,7 @@ namespace Opm
std::vector<MAT_SIZE_T> piv(num_basis); std::vector<MAT_SIZE_T> piv(num_basis);
MAT_SIZE_T ldb = num_basis; MAT_SIZE_T ldb = num_basis;
MAT_SIZE_T info = 0; MAT_SIZE_T info = 0;
orig_jac_ = jac_;
dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info); dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info);
if (info != 0) { if (info != 0) {
// Print the local matrix and rhs. // Print the local matrix and rhs.
@ -624,7 +626,7 @@ namespace Opm
<< " with A = \n"; << " with A = \n";
for (int row = 0; row < n; ++row) { for (int row = 0; row < n; ++row) {
for (int col = 0; col < n; ++col) { for (int col = 0; col < n; ++col) {
std::cerr << " " << jac_[row + n*col]; std::cerr << " " << orig_jac_[row + n*col];
} }
std::cerr << '\n'; std::cerr << '\n';
} }

1
opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
View File

@ -91,6 +91,7 @@ namespace Opm
double* tof_coeff_; double* tof_coeff_;
std::vector<double> rhs_; // single-cell right-hand-side std::vector<double> rhs_; // single-cell right-hand-side
std::vector<double> jac_; // single-cell jacobian std::vector<double> jac_; // single-cell jacobian
std::vector<double> orig_jac_; // single-cell jacobian (copy)
// Below: storage for quantities needed by solveSingleCell(). // Below: storage for quantities needed by solveSingleCell().
std::vector<double> coord_; std::vector<double> coord_;
std::vector<double> basis_; std::vector<double> basis_;

4
opm/core/utility/VelocityInterpolation.cpp
View File

@ -178,6 +178,7 @@ namespace Opm
// for each corner. // for each corner.
const int dim = grid_.dimensions; const int dim = grid_.dimensions;
std::vector<double> N(dim*dim); // Normals matrix. Fortran ordering! std::vector<double> N(dim*dim); // Normals matrix. Fortran ordering!
std::vector<double> orig_N(dim*dim); // Normals matrix. Fortran ordering!
std::vector<double> f(dim); // Flux vector. std::vector<double> f(dim); // Flux vector.
std::vector<MAT_SIZE_T> piv(dim); // For LAPACK solve std::vector<MAT_SIZE_T> piv(dim); // For LAPACK solve
const int num_cells = grid_.number_of_cells; const int num_cells = grid_.number_of_cells;
@ -203,6 +204,7 @@ namespace Opm
MAT_SIZE_T lda = n; MAT_SIZE_T lda = n;
MAT_SIZE_T ldb = n; MAT_SIZE_T ldb = n;
MAT_SIZE_T info = 0; MAT_SIZE_T info = 0;
orig_N = N;
dgesv_(&n, &nrhs, &N[0], &lda, &piv[0], &f[0], &ldb, &info); dgesv_(&n, &nrhs, &N[0], &lda, &piv[0], &f[0], &ldb, &info);
if (info != 0) { if (info != 0) {
// Print the local matrix and rhs. // Print the local matrix and rhs.
@ -210,7 +212,7 @@ namespace Opm
<< " with N = \n"; << " with N = \n";
for (int row = 0; row < n; ++row) { for (int row = 0; row < n; ++row) {
for (int col = 0; col < n; ++col) { for (int col = 0; col < n; ++col) {
std::cerr << " " << N[row + n*col]; std::cerr << " " << orig_N[row + n*col];
} }
std::cerr << '\n'; std::cerr << '\n';
} }