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

merge

Browse Source
This commit is contained in:
Kjetil Olsen Lye
2012-03-26 11:17:36 +02:00
parent dd2d474643 5c2cc54ecc
commit 3f63435100
9 changed files with 234 additions and 149 deletions
Download Patch File Download Diff File Expand all files Collapse all files

214
opm/core/eclipse/EclipseGridParser.cpp
View File

@@ -170,6 +170,98 @@ namespace {
}
}
inline std::string upcase(const std::string& s)
{
std::string us(s);
// Getting the character type facet for toupper().
// We use the classic (i.e. C) locale.
const std::ctype<char>& ct = std::use_facet< std::ctype<char> >(std::locale::classic());
for (int i = 0; i < int(s.size()); ++i) {
us[i] = ct.toupper(s[i]);
}
return us;
}
inline std::string readKeyword(std::istream& is)
{
std::string keyword_candidate;
while (!is.eof()) {
is >> keyword_candidate;
if(keyword_candidate.find("--") == 0) {
is >> ignoreLine; // This line is a comment
} else {
return upcase(keyword_candidate);
}
}
return "CONTINUE"; // Last line in included file is a comment
}
inline bool readKeywordNew(std::istream& is, std::string& keyword)
{
char buf[9];
int i, j;
char c;
/* Clear buf */
for (i=0; i<9; ++i) {
buf[i] = '\0';
}
/* Read first character and check if it is uppercase*/
//buf[0] = fgetc(fp);
is.get(buf[0]);
if ( !isupper( buf[0] ) ) {
is.unget();
return false; /* NOT VALID CHARACTER */
}
/* Scan as much as possible possible keyword, 8 characters long */
i = 1;
is.get(c);
while ( (is.good()) &&
(c != EOF ) &&
(!isblank(c) ) &&
(isupper(c) || isdigit(c)) &&
(c != '\n' ) &&
(c != '/' ) &&
(i < 8 )) {
buf[i++] = c;
is.get(c);
}
/* Skip rest of line */
if (c != '\n'){
is.get(c);
while ( (is.good()) &&
(c != EOF ) &&
(c != '\n' )) {
is.get(c);
}
}
if(c == '\n') {
is.unget();
}
/* Find first non-uppercase or non-digit character */
for (i=0; i<8; ++i) {
if ( !(isupper(buf[i]) || isdigit(buf[i])) ) {
break;
}
}
/* Check if remaining characters are blank */
for (j = i; j<8; ++j) {
if(!isspace(buf[j]) && buf[j] != '\0') {
return false; /* CHARACTER AFTER SPACE OR INVALID CHARACTER */
}
buf[j] = '\0';
}
keyword = std::string(buf);
std::string::size_type end = keyword.find_last_of('\0');
if(end != keyword.npos)
keyword = keyword.substr(0, end+1);
return true;
}
} // anon namespace
@@ -228,7 +320,7 @@ void EclipseGridParser::readImpl(istream& is)
// Make temporary maps that will at the end be swapped with the
// member maps
// NOTE: Above is no longer true, for easier implementation of
// NOTE: Above is no longer true, for easier implementation of
// the INCLUDE keyword. We lose the strong exception guarantee,
// though (of course retaining the basic guarantee).
map<string, vector<int> >& intmap = integer_field_map_;
@@ -236,71 +328,71 @@ void EclipseGridParser::readImpl(istream& is)
map<string, tr1::shared_ptr<SpecialBase> >& specialmap = special_field_map_;
// Actually read the data
is >> ignoreWhitespace;
while (!is.eof()) {
string keyword = readKeyword(is);
#ifdef VERBOSE
cout << "Keyword found: " << keyword << endl;
#endif
FieldType type = classifyKeyword(keyword);
switch (type) {
case Integer:
readVectorData(is, intmap[keyword]);
break;
case FloatingPoint:
readVectorData(is, floatmap[keyword]);
break;
case SpecialField: {
map<string, std::tr1::shared_ptr<SpecialBase> >::iterator pos =
specialmap.find(keyword);
if (pos == specialmap.end()) {
std::tr1::shared_ptr<SpecialBase> sb_ptr =
createSpecialField(is, keyword);
std::string keyword;
while (is.good()) {
is >> ignoreWhitespace;
bool ok = readKeywordNew(is, keyword);
if (ok) {
//#ifdef VERBOSE
cout << "Keyword found: " << keyword << endl;
//#endif
FieldType type = classifyKeyword(keyword);
switch (type) {
case Integer:
readVectorData(is, intmap[keyword]);
break;
case FloatingPoint:
readVectorData(is, floatmap[keyword]);
break;
case SpecialField: {
std::tr1::shared_ptr<SpecialBase> sb_ptr = createSpecialField(is, keyword);
if (sb_ptr) {
specialmap[keyword] = sb_ptr;
} else {
THROW("Could not create field " << keyword);
}
} else {
pos->second->read(is);
break;
}
break;
}
case IgnoreWithData: {
ignored_fields_.insert(keyword);
is >> ignoreSlashLine;
#ifdef VERBOSE
cout << "(ignored)" << endl;
#endif
break;
}
case IgnoreNoData: {
ignored_fields_.insert(keyword);
case IgnoreWithData: {
ignored_fields_.insert(keyword);
//is >> ignoreSlashLine;
//#ifdef VERBOSE
// cout << "(ignored)" << endl;
//#endif
break;
}
case IgnoreNoData: {
ignored_fields_.insert(keyword);
//is >> ignoreLine;
//#ifdef VERBOSE
// cout << "(ignored)" << endl;
//#endif
break;
}
case Include: {
string include_filename = readString(is);
if (!directory_.empty()) {
include_filename = directory_ + '/' + include_filename;
}
ifstream include_is(include_filename.c_str());
if (!include_is) {
THROW("Unable to open INCLUDEd file " << include_filename);
}
readImpl(include_is);
// is >> ignoreSlashLine;
break;
}
case Unknown:
default:
ignored_fields_.insert(keyword);
cout << "*** Warning: keyword " << keyword << " is unknown." << endl;
//is >> ignoreSlashLine;
//throw exception();
}
} else {
// if (!ok)
is >> ignoreLine;
#ifdef VERBOSE
cout << "(ignored)" << endl;
#endif
break;
}
case Include: {
string include_filename = readString(is);
if (!directory_.empty()) {
include_filename = directory_ + '/' + include_filename;
}
ifstream include_is(include_filename.c_str());
if (!include_is) {
THROW("Unable to open INCLUDEd file " << include_filename);
}
readImpl(include_is);
is >> ignoreSlashLine;
break;
}
case Unknown:
default:
cerr << "Keyword " << keyword << " not recognized." << endl;
throw exception();
}
is >> ignoreWhitespace;
}
#define VERBOSE_LIST_FIELDS 0
@@ -356,7 +448,7 @@ void EclipseGridParser::convertToSI()
unit = 1.0;
do_convert = false; // Dimensionless keywords...
} else if (key == "PRESSURE") {
unit = units_.pressure;
unit = units_.pressure;
} else {
THROW("Units for field " << key << " not specified. Cannon convert to SI.");
}
@@ -452,7 +544,7 @@ const std::vector<int>& EclipseGridParser::getIntegerValue(const std::string& ke
map<string, vector<int> >::const_iterator it
= integer_field_map_.find(keyword);
if (it == integer_field_map_.end()) {
return empty_integer_field_;
THROW("No such field: " << keyword);
} else {
return it->second;
}
@@ -465,7 +557,7 @@ const std::vector<double>& EclipseGridParser::getFloatingPointValue(const std::s
map<string, vector<double> >::const_iterator it
= floating_field_map_.find(keyword);
if (it == floating_field_map_.end()) {
return empty_floating_field_;
THROW("No such field: " << keyword);
} else {
return it->second;
}

2
opm/core/eclipse/EclipseGridParser.hpp
View File

@@ -186,8 +186,6 @@ private:
std::map<std::string, std::vector<double> > floating_field_map_;
std::map<std::string, std::tr1::shared_ptr<SpecialBase> > special_field_map_;
std::set<std::string> ignored_fields_;
std::vector<int> empty_integer_field_;
std::vector<double> empty_floating_field_;
std::tr1::shared_ptr<SpecialBase> empty_special_field_;
EclipseUnits units_;
};

26
opm/core/eclipse/EclipseGridParserHelpers.hpp
View File

@@ -78,32 +78,6 @@ namespace
return is;
}
inline std::string upcase(const std::string& s)
{
std::string us(s);
// Getting the character type facet for toupper().
// We use the classic (i.e. C) locale.
const std::ctype<char>& ct = std::use_facet< std::ctype<char> >(std::locale::classic());
for (int i = 0; i < int(s.size()); ++i) {
us[i] = ct.toupper(s[i]);
}
return us;
}
inline std::string readKeyword(std::istream& is)
{
std::string keyword_candidate;
while (!is.eof()) {
is >> keyword_candidate;
if(keyword_candidate.find("--") == 0) {
is >> ignoreLine; // This line is a comment
} else {
return upcase(keyword_candidate);
}
}
return "CONTINUE"; // Last line in included file is a comment
}
inline std::string readString(std::istream& is)
{
std::string string_candidate;

10
opm/core/linalg/blas_lapack.h
View File

@@ -94,6 +94,16 @@ void dgbsv_(const MAT_SIZE_T *n ,
const MAT_SIZE_T *ldb ,
MAT_SIZE_T *info);
/* B <- A \ B, general solver */
void dgesv_(const MAT_SIZE_T *n,
const MAT_SIZE_T *nrhs ,
double *A ,
const MAT_SIZE_T *lda ,
MAT_SIZE_T *piv ,
double *B ,
const MAT_SIZE_T *ldb ,
MAT_SIZE_T *info);
/* A <- chol(A) */
void dpotrf_(const char *uplo, const MAT_SIZE_T *n,
double *A , const MAT_SIZE_T *lda,

2
opm/core/pressure/HybridPressureSolver.hpp
View File

@@ -129,7 +129,7 @@ public:
const std::vector<double>& total_mobilities,
const std::vector<double>& omegas,
const std::vector<FlowBCTypes>& bctypes,
const std::vector<double> bcvalues)
const std::vector<double>& bcvalues)
{
if (state_ == Uninitialized) {
throw std::runtime_error("Error in HybridPressureSolver::assemble(): You must call init() prior to calling assemble().");

2
opm/core/pressure/TPFAPressureSolver.hpp
View File

@@ -115,7 +115,7 @@ public:
const std::vector<double>& total_mobilities,
const std::vector<double>& omegas,
const std::vector<FlowBCTypes>& bctypes,
const std::vector<double> bcvalues)
const std::vector<double>& bcvalues)
{
if (state_ == Uninitialized) {
throw std::runtime_error("Error in TPFAPressureSolver::assemble(): You must call init() prior to calling assemble().");

32
opm/core/utility/parameters/ParameterGroup.cpp
View File

@@ -119,19 +119,25 @@ namespace Opm {
int lineno = 0;
while (samcode_readline(is, parameter)) {
++lineno;
int fpos = parameter.find(ID_delimiter_assignment);
if (fpos == int(std::string::npos)) {
std::cerr << "WARNING: No '" << ID_delimiter_assignment << "' found on line " << lineno << ".\n";
}
int pos = fpos + ID_delimiter_assignment.size();
int spos = parameter.find(ID_delimiter_assignment, pos);
if (spos == int(std::string::npos)) {
std::string name = parameter.substr(0, fpos);
std::string value = parameter.substr(pos, spos);
this->insertParameter(name, value);
} else {
std::cerr << "WARNING: To many '" << ID_delimiter_assignment << "' found on line " << lineno << ".\n";
}
int commentpos = parameter.find(ID_comment);
if (commentpos != 0) {
if (commentpos != int(std::string::npos)) {
parameter = parameter.substr(0, commentpos);
}
int fpos = parameter.find(ID_delimiter_assignment);
if (fpos == int(std::string::npos)) {
std::cerr << "WARNING: No '" << ID_delimiter_assignment << "' found on line " << lineno << ".\n";
}
int pos = fpos + ID_delimiter_assignment.size();
int spos = parameter.find(ID_delimiter_assignment, pos);
if (spos == int(std::string::npos)) {
std::string name = parameter.substr(0, fpos);
std::string value = parameter.substr(pos, spos);
this->insertParameter(name, value);
} else {
std::cerr << "WARNING: To many '" << ID_delimiter_assignment << "' found on line " << lineno << ".\n";
}
}
}
#ifdef MATLAB_MEX_FILE
fclose(is);

1
opm/core/utility/parameters/ParameterStrings.hpp
View File

@@ -63,6 +63,7 @@ namespace Opm {
const std::string ID_path_root = "";
const std::string ID_delimiter_path = "/";
const std::string ID_comment = "//";
const std::string ID_delimiter_assignment = "=";
} // namespace parameter
} // namespace Opm

94
tests/test_lapack.cpp
View File

@@ -61,77 +61,81 @@ int main()
} else {
std::cerr << "Something went wrong in dgtsv_()\n";
}
std::cout << std::endl;
//test of dgbsv_
int ldb = N;
int lda = N;
std::vector<int> ipiv(N, 0);
std::vector<double> AA(N*N, 0.);
std::vector<double> BB(N, 0.);
for (int i = 0; i < N; ++i) {
AA[i + i*N] = 10.;
if (i > 0) {
AA[i + (i - 1)*N] = i;
AA[i - 1 + i*N] = N - i;
}
BB[i] = i;
}
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
std::cout << " " << AA[i + j*N];
}
std::cout << " " << std::endl;
}
std::cout << std::endl;
int kl = 1;
int ku = 1;
int nrowAB = 2*kl + ku + 1;
int ldb = N;
int ldab = nrowAB;
int ipiv;
std::vector<double> AB(nrowAB*N, 0.);
std::vector<double> BB(N, 0.);
double ABu[N] = { 0., 2.1, -1.0, 1.9, 8.0};
double ABd[N] = { 3.0, 2.3, -5.0, -0.9, 7.1 };
double ABl[N] = { 3.4, 3.6, 7.0, - 6.0, 0. };
for (int i = 0; i<N; ++i) {
AB[kl + i*nrowAB] = ABu[i];
AB[kl + 1 + i*nrowAB] = ABd[i];
AB[kl + 2 + i*nrowAB] = ABl[i];
BandMatrixCoeff bmc(N, ku, kl);
// Rewrite AA matrix in band format AB
for (int i = 0; i < N; ++i) {
for (int j = -1; j < 2; ++j) {
if (i + j > -1 && i + j < N)
AB[bmc(i, i + j)] = AA[i + N*(i + j)];
}
}
BB[0] = 2.7;
BB[1] = -0.5;
BB[2] = 2.6;
BB[3] = 0.6;
BB[4] = 2.7;
for (int i = 0; i < nrowAB; ++i) {
for (int j = 0; j < N; ++j) {
std::cout << " " << AB[i + j*nrowAB];
}
std::cout << " " << std::endl;
}
std::cout << std::endl;
dgbsv_(&N, &kl, &ku, &nrhs, &AB[0], &ldab, &ipiv, &BB[0], &ldb, &info);
dgesv_(&N ,&nrhs, &AA[0], &lda, &ipiv[0], &BB[0], &ldb, &info);
if (info == 0) {
for (int i = 0; i < N; ++i) {
std::cout << BB[i] << ' ';
}
std::cout << std::endl;
} else {
std::cerr << "Something went wrong in dgbsv_()\n";
std::cerr << "Something went wrong in debsv_()\n";
}
int NN = 3;
kl = 1;
ku = 1;
nrowAB = 2*kl + ku + 1;
ldb = NN;
ldab = nrowAB;
AB.assign(NN*nrowAB, 0.);
BB.assign(NN, 0.);
BandMatrixCoeff bmc(NN, ku, kl);
AB[bmc(0, 0)] = 1.;
AB[bmc(1, 0)] = 1.;
AB[bmc(0, 1)] = 1.;
AB[bmc(1, 1)] = 2.;
AB[bmc(2, 1)] = 2.;
AB[bmc(1, 2)] = 1.;
AB[bmc(2, 2)] = 4.;
BB[0] = 3.;
BB[1] = 8.;
BB[2] = 16.;
for (int i = 0; i < N; ++i) {
BB[i] = i;
}
// std::vector<double>::iterator it;
// std::cout << "myvector contains:";
// for ( it=AB.begin() ; it < AB.end(); it++ ) {
// std::cout << " " << *it;
// }
// std::cout << std::endl;
dgbsv_(&N, &kl, &ku, &nrhs, &AB[0], &ldab, &ipiv[0], &BB[0], &ldb, &info);
dgbsv_(&NN ,&kl, &ku, &nrhs, &AB[0], &ldab, &ipiv, &BB[0], &ldb, &info);
if (info == 0) {
for (int i = 0; i < NN; ++i) {
for (int i = 0; i < N; ++i) {
std::cout << BB[i] << ' ';
}
std::cout << std::endl;
} else {
std::cerr << "Something went wrong in dgbsv_()\n";
std::cerr << "Something went wrong in debsv_()\n";
}
}