From 7ec7651ea35da28323a09a6bbac9077c22c141b5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 5 Jun 2012 11:54:03 +0200 Subject: [PATCH 01/23] Work in progress on adding some SCHEDULE support. --- opm/core/eclipse/EclipseGridParser.cpp | 136 ++++++++++++++++++------- opm/core/eclipse/EclipseGridParser.hpp | 27 ++++- 2 files changed, 122 insertions(+), 41 deletions(-) diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index 29d8bb05..e7e7f9cc 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -147,6 +147,7 @@ namespace { enum FieldType { Integer, FloatingPoint, + Timestepping, SpecialField, IgnoreWithData, IgnoreNoData, @@ -161,6 +162,8 @@ namespace { return Integer; } else if (count(floating_fields, floating_fields + num_floating_fields, keyword)) { return FloatingPoint; + } else if (keyword == "TSTEP" || keyword == "DATES") { + return Timestepping; } else if (count(special_fields, special_fields + num_special_fields, keyword)) { return SpecialField; } else if (count(ignore_with_data, ignore_with_data + num_ignore_with_data, keyword)) { @@ -276,6 +279,8 @@ namespace { //--------------------------------------------------------------------------- EclipseGridParser::EclipseGridParser() //--------------------------------------------------------------------------- + : current_reading_mode_(NonTimesteps), + current_epoch_(0) { } @@ -284,6 +289,8 @@ EclipseGridParser::EclipseGridParser() //--------------------------------------------------------------------------- EclipseGridParser::EclipseGridParser(const string& filename, bool convert_to_SI) //--------------------------------------------------------------------------- + : current_reading_mode_(NonTimesteps), + current_epoch_(0) { // Store directory of filename boost::filesystem::path p(filename); @@ -304,13 +311,39 @@ void EclipseGridParser::read(istream& is, bool convert_to_SI) { integer_field_map_.clear(); floating_field_map_.clear(); - special_field_map_.clear(); + special_field_by_epoch_.clear(); + special_field_by_epoch_.push_back(SpecialMap()); readImpl(is); + + current_epoch_ = 0; + computeUnits(); if (convert_to_SI) { convertToSI(); } + +#define VERBOSE_LIST_FIELDS 1 +#if VERBOSE_LIST_FIELDS + std::cout << "\nInteger fields:\n"; + for (std::map >::iterator + i = integer_field_map_.begin(); i != integer_field_map_.end(); ++i) + std::cout << '\t' << i->first << '\n'; + + std::cout << "\nFloat fields:\n"; + for (std::map >::iterator + i = floating_field_map_.begin(); i != floating_field_map_.end(); ++i) + std::cout << '\t' << i->first << '\n'; + + std::cout << "\nSpecial fields:\n"; + for (int epoch = 0; epoch < numberOfEpochs(); ++epoch) { + std::cout << "Epoch " << epoch << '\n'; + const SpecialMap& sm = special_field_by_epoch_[epoch]; + for (SpecialMap::const_iterator i = sm.begin(); i != sm.end(); ++i) { + std::cout << '\t' << i->first << '\n'; + } + } +#endif } //--------------------------------------------------------------------------- @@ -329,7 +362,6 @@ void EclipseGridParser::readImpl(istream& is) // though (of course retaining the basic guarantee). map >& intmap = integer_field_map_; map >& floatmap = floating_field_map_; - map >& specialmap = special_field_map_; // Actually read the data std::string keyword; @@ -341,6 +373,7 @@ void EclipseGridParser::readImpl(istream& is) cout << "Keyword found: " << keyword << endl; //#endif FieldType type = classifyKeyword(keyword); + // std::cout << "Classification: " << type << std::endl; switch (type) { case Integer: readVectorData(is, intmap[keyword]); @@ -348,14 +381,42 @@ void EclipseGridParser::readImpl(istream& is) case FloatingPoint: readVectorData(is, floatmap[keyword]); break; - case SpecialField: { - std::tr1::shared_ptr sb_ptr = createSpecialField(is, keyword); - if (sb_ptr) { - specialmap[keyword] = sb_ptr; - } else { - THROW("Could not create field " << keyword); + case Timestepping: + if (current_reading_mode_ == NonTimesteps) { + current_reading_mode_ = Timesteps; + } + // Append to current epoch's TSTEP. + // Update current_reading_date_? + // Currently do neither... + { // A scope to isolate sb_ptr. + SpecialFieldPtr sb_ptr = createSpecialField(is, keyword); + if (sb_ptr) { + special_field_by_epoch_[current_epoch_][keyword] = sb_ptr; + } else { + THROW("Could not create field " << keyword); + } } break; + case SpecialField: { + if (current_reading_mode_ == Timesteps) { + // We have been reading timesteps, but have + // now encountered something else. + // That means we are in a new epoch. + current_reading_mode_ = NonTimesteps; + // New epoch starts out as a copy of old epoch. + special_field_by_epoch_.push_back(special_field_by_epoch_.back()); + ++current_epoch_; + ASSERT(int(special_field_by_epoch_.size()) == current_epoch_ + 1); + } + { // A scope to isolate sb_ptr. + SpecialFieldPtr sb_ptr = createSpecialField(is, keyword); + if (sb_ptr) { + special_field_by_epoch_[current_epoch_][keyword] = sb_ptr; + } else { + THROW("Could not create field " << keyword); + } + break; + } } case IgnoreWithData: { ignored_fields_.insert(keyword); @@ -398,24 +459,6 @@ void EclipseGridParser::readImpl(istream& is) is >> ignoreLine; } } - -#define VERBOSE_LIST_FIELDS 0 -#if VERBOSE_LIST_FIELDS - std::cout << "\nInteger fields:\n"; - for (std::map >::iterator - i = intmap.begin(); i != intmap.end(); ++i) - std::cout << '\t' << i->first << '\n'; - - std::cout << "\nFloat fields:\n"; - for (std::map >::iterator - i = floatmap.begin(); i != floatmap.end(); ++i) - std::cout << '\t' << i->first << '\n'; - - std::cout << "\nSpecial fields:\n"; - for (std::map >::iterator - i = specialmap.begin(); i != specialmap.end(); ++i) - std::cout << '\t' << i->first << '\n'; -#endif } @@ -425,9 +468,12 @@ void EclipseGridParser::convertToSI() //--------------------------------------------------------------------------- { // Convert all special fields. - typedef std::map >::iterator SpecialIt; - for (SpecialIt i = special_field_map_.begin(); i != special_field_map_.end(); ++i) { - i->second->convertToSI(units_); + typedef SpecialMap::iterator SpecialIt; + for (int epoch = 0; epoch < numberOfEpochs(); ++epoch) { + SpecialMap& sm = special_field_by_epoch_[epoch]; + for (SpecialIt i = sm.begin(); i != sm.end(); ++i) { + i->second->convertToSI(units_); + } } // Convert all floating point fields. @@ -479,7 +525,7 @@ bool EclipseGridParser::hasField(const string& keyword) const { string ukey = upcase(keyword); return integer_field_map_.count(ukey) || floating_field_map_.count(ukey) || - special_field_map_.count(ukey) || ignored_fields_.count(ukey); + special_field_by_epoch_[current_epoch_].count(ukey) || ignored_fields_.count(ukey); } @@ -505,7 +551,7 @@ vector EclipseGridParser::fieldNames() const vector names; names.reserve(integer_field_map_.size() + floating_field_map_.size() + - special_field_map_.size() + + special_field_by_epoch_[current_epoch_].size() + ignored_fields_.size()); { map >::const_iterator it = integer_field_map_.begin(); @@ -520,8 +566,8 @@ vector EclipseGridParser::fieldNames() const } } { - map >::const_iterator it = special_field_map_.begin(); - for (; it != special_field_map_.end(); ++it) { + SpecialMap::const_iterator it = special_field_by_epoch_[current_epoch_].begin(); + for (; it != special_field_by_epoch_[current_epoch_].end(); ++it) { names.push_back(it->first); } } @@ -534,6 +580,24 @@ vector EclipseGridParser::fieldNames() const return names; } + +//--------------------------------------------------------------------------- +int EclipseGridParser::numberOfEpochs() const +//--------------------------------------------------------------------------- +{ + return special_field_by_epoch_.size(); +} + + +//--------------------------------------------------------------------------- +void EclipseGridParser::setCurrentEpoch(int epoch) +//--------------------------------------------------------------------------- +{ + ASSERT(epoch >= 0 && epoch < numberOfEpochs()); + current_epoch_ = epoch; +} + + //--------------------------------------------------------------------------- const std::vector& EclipseGridParser::getIntegerValue(const std::string& keyword) const //--------------------------------------------------------------------------- @@ -574,8 +638,8 @@ const std::vector& EclipseGridParser::getFloatingPointValue(const std::s const std::tr1::shared_ptr EclipseGridParser::getSpecialValue(const std::string& keyword) const //--------------------------------------------------------------------------- { - map >::const_iterator it = special_field_map_.find(keyword); - if (it == special_field_map_.end()) { + SpecialMap::const_iterator it = special_field_by_epoch_[current_epoch_].find(keyword); + if (it == special_field_by_epoch_[current_epoch_].end()) { THROW("No such field: " << keyword); } else { return it->second; @@ -617,7 +681,7 @@ void EclipseGridParser::setSpecialField(const std::string& keyword, std::tr1::shared_ptr field) //--------------------------------------------------------------------------- { - special_field_map_[keyword] = field; + special_field_by_epoch_[current_epoch_][keyword] = field; } //--------------------------------------------------------------------------- diff --git a/opm/core/eclipse/EclipseGridParser.hpp b/opm/core/eclipse/EclipseGridParser.hpp index 5662ae87..4ddac173 100644 --- a/opm/core/eclipse/EclipseGridParser.hpp +++ b/opm/core/eclipse/EclipseGridParser.hpp @@ -94,6 +94,14 @@ public: /// The keywords/fields found in the file. std::vector fieldNames() const; + /// Returns the number of distinct epochs found in the deck SCHEDULE. + int numberOfEpochs() const; + + /// Sets the current epoch. + /// Valid arguments are in [0, ..., numberOfEpochs() - 1]. + /// After reading, current epoch always starts at 0. + void setCurrentEpoch(int epoch); + /// Returns a reference to a vector containing the values /// corresponding to the given integer keyword. const std::vector& getIntegerValue(const std::string& keyword) const; @@ -102,10 +110,12 @@ public: /// corresponding to the given floating-point keyword. const std::vector& getFloatingPointValue(const std::string& keyword) const; + typedef std::tr1::shared_ptr SpecialFieldPtr; + /// Returns a reference to a vector containing pointers to the values /// corresponding to the given keyword when the values are not only integers /// or floats. - const std::tr1::shared_ptr getSpecialValue(const std::string& keyword) const; + const SpecialFieldPtr getSpecialValue(const std::string& keyword) const; // This macro implements support for a special field keyword. It requires that a subclass // of SpecialBase exists, that has the same name as the keyword. @@ -170,7 +180,7 @@ public: void setFloatingPointField(const std::string& keyword, const std::vector& field); /// Sets a special field to have a particular value. - void setSpecialField(const std::string& keyword, std::tr1::shared_ptr field); + void setSpecialField(const std::string& keyword, SpecialFieldPtr field); /// Compute the units used by the deck, depending on the presence /// of keywords such as METRIC, FIELD etc. It is an error to call @@ -181,17 +191,24 @@ public: const EclipseUnits& units() const; private: - std::tr1::shared_ptr createSpecialField(std::istream& is, const std::string& fieldname); + SpecialFieldPtr createSpecialField(std::istream& is, const std::string& fieldname); void readImpl(std::istream& is); std::string directory_; std::map > integer_field_map_; std::map > floating_field_map_; - std::map > special_field_map_; + // std::map special_field_map_; std::set ignored_fields_; - std::tr1::shared_ptr empty_special_field_; EclipseUnits units_; + + // For SCHEDULE handling. + enum ReadingMode { NonTimesteps, Timesteps }; + ReadingMode current_reading_mode_; + boost::gregorian::date current_reading_date_; + int current_epoch_; + typedef std::map SpecialMap; + std::vector special_field_by_epoch_; }; From 9a2adee21afaeebf703f2b3c5d1585b54e72a6b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 5 Jun 2012 13:02:47 +0200 Subject: [PATCH 02/23] Allow specifying well rates to be zero. --- opm/core/WellsManager.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/opm/core/WellsManager.cpp b/opm/core/WellsManager.cpp index 6d0d83dd..8aa5d644 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/WellsManager.cpp @@ -435,13 +435,13 @@ namespace Opm // Add all controls that are present in well. int ok = 1; int control_pos[5] = { -1, -1, -1, -1, -1 }; - if (ok && wci_line.surface_flow_max_rate_ > 0.0) { + if (ok && wci_line.surface_flow_max_rate_ >= 0.0) { control_pos[InjectionControl::RATE] = w_->ctrls[wix]->num; const double distr[3] = { 1.0, 1.0, 1.0 }; ok = append_well_controls(SURFACE_RATE, wci_line.surface_flow_max_rate_, distr, wix, w_); } - if (ok && wci_line.reservoir_flow_max_rate_ > 0.0) { + if (ok && wci_line.reservoir_flow_max_rate_ >= 0.0) { control_pos[InjectionControl::RESV] = w_->ctrls[wix]->num; const double distr[3] = { 1.0, 1.0, 1.0 }; ok = append_well_controls(RESERVOIR_RATE, wci_line.reservoir_flow_max_rate_, @@ -515,7 +515,7 @@ namespace Opm // Add all controls that are present in well. int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 }; int ok = 1; - if (ok && wcp_line.oil_max_rate_ > 0.0) { + if (ok && wcp_line.oil_max_rate_ >= 0.0) { if (!pu.phase_used[BlackoilPhases::Liquid]) { THROW("Oil phase not active and ORAT control specified."); } @@ -525,7 +525,7 @@ namespace Opm ok = append_well_controls(SURFACE_RATE, -wcp_line.oil_max_rate_, distr, wix, w_); } - if (ok && wcp_line.water_max_rate_ > 0.0) { + if (ok && wcp_line.water_max_rate_ >= 0.0) { if (!pu.phase_used[BlackoilPhases::Aqua]) { THROW("Water phase not active and WRAT control specified."); } @@ -535,7 +535,7 @@ namespace Opm ok = append_well_controls(SURFACE_RATE, -wcp_line.water_max_rate_, distr, wix, w_); } - if (ok && wcp_line.gas_max_rate_ > 0.0) { + if (ok && wcp_line.gas_max_rate_ >= 0.0) { if (!pu.phase_used[BlackoilPhases::Vapour]) { THROW("Gas phase not active and GRAT control specified."); } @@ -545,7 +545,7 @@ namespace Opm ok = append_well_controls(SURFACE_RATE, -wcp_line.gas_max_rate_, distr, wix, w_); } - if (ok && wcp_line.liquid_max_rate_ > 0.0) { + if (ok && wcp_line.liquid_max_rate_ >= 0.0) { if (!pu.phase_used[BlackoilPhases::Aqua]) { THROW("Water phase not active and LRAT control specified."); } @@ -559,7 +559,7 @@ namespace Opm ok = append_well_controls(SURFACE_RATE, -wcp_line.liquid_max_rate_, distr, wix, w_); } - if (ok && wcp_line.reservoir_flow_max_rate_ > 0.0) { + if (ok && wcp_line.reservoir_flow_max_rate_ >= 0.0) { control_pos[ProductionControl::RESV] = w_->ctrls[wix]->num; double distr[3] = { 1.0, 1.0, 1.0 }; ok = append_well_controls(RESERVOIR_RATE, -wcp_line.reservoir_flow_max_rate_, From f1bc8ca160fdddd371c170fc376f0e651b55ba4f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 5 Jun 2012 13:04:57 +0200 Subject: [PATCH 03/23] Work in progress on SCHEDULE. Now handles TSTEP properly. --- opm/core/eclipse/EclipseGridParser.cpp | 47 ++++++++++++++------------ opm/core/eclipse/EclipseGridParser.hpp | 2 +- 2 files changed, 26 insertions(+), 23 deletions(-) diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index e7e7f9cc..0780faaf 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -279,7 +279,7 @@ namespace { //--------------------------------------------------------------------------- EclipseGridParser::EclipseGridParser() //--------------------------------------------------------------------------- - : current_reading_mode_(NonTimesteps), + : current_reading_mode_(Regular), current_epoch_(0) { } @@ -289,7 +289,7 @@ EclipseGridParser::EclipseGridParser() //--------------------------------------------------------------------------- EclipseGridParser::EclipseGridParser(const string& filename, bool convert_to_SI) //--------------------------------------------------------------------------- - : current_reading_mode_(NonTimesteps), + : current_reading_mode_(Regular), current_epoch_(0) { // Store directory of filename @@ -375,48 +375,51 @@ void EclipseGridParser::readImpl(istream& is) FieldType type = classifyKeyword(keyword); // std::cout << "Classification: " << type << std::endl; switch (type) { - case Integer: + case Integer: { readVectorData(is, intmap[keyword]); break; - case FloatingPoint: + } + case FloatingPoint: { readVectorData(is, floatmap[keyword]); break; - case Timestepping: - if (current_reading_mode_ == NonTimesteps) { + } + case Timestepping: { + if (current_reading_mode_ == Regular) { current_reading_mode_ = Timesteps; } - // Append to current epoch's TSTEP. - // Update current_reading_date_? - // Currently do neither... - { // A scope to isolate sb_ptr. + SpecialMap& sm = special_field_by_epoch_[current_epoch_]; + // Append to current epoch's TSTEP, if it exists. + SpecialMap::iterator it = sm.find("TSTEP"); + if (it != sm.end()) { + it->second->read(is); // This will append to the TSTEP object. + } else { + // There is no existing TSTEP for this epoch, create it. SpecialFieldPtr sb_ptr = createSpecialField(is, keyword); if (sb_ptr) { - special_field_by_epoch_[current_epoch_][keyword] = sb_ptr; + sm[keyword] = sb_ptr; } else { THROW("Could not create field " << keyword); } } break; + } case SpecialField: { if (current_reading_mode_ == Timesteps) { // We have been reading timesteps, but have // now encountered something else. // That means we are in a new epoch. - current_reading_mode_ = NonTimesteps; - // New epoch starts out as a copy of old epoch. - special_field_by_epoch_.push_back(special_field_by_epoch_.back()); + current_reading_mode_ = Regular; + special_field_by_epoch_.push_back(SpecialMap()); ++current_epoch_; ASSERT(int(special_field_by_epoch_.size()) == current_epoch_ + 1); } - { // A scope to isolate sb_ptr. - SpecialFieldPtr sb_ptr = createSpecialField(is, keyword); - if (sb_ptr) { - special_field_by_epoch_[current_epoch_][keyword] = sb_ptr; - } else { - THROW("Could not create field " << keyword); - } - break; + SpecialFieldPtr sb_ptr = createSpecialField(is, keyword); + if (sb_ptr) { + special_field_by_epoch_[current_epoch_][keyword] = sb_ptr; + } else { + THROW("Could not create field " << keyword); } + break; } case IgnoreWithData: { ignored_fields_.insert(keyword); diff --git a/opm/core/eclipse/EclipseGridParser.hpp b/opm/core/eclipse/EclipseGridParser.hpp index 4ddac173..cf57906b 100644 --- a/opm/core/eclipse/EclipseGridParser.hpp +++ b/opm/core/eclipse/EclipseGridParser.hpp @@ -203,7 +203,7 @@ private: EclipseUnits units_; // For SCHEDULE handling. - enum ReadingMode { NonTimesteps, Timesteps }; + enum ReadingMode { Regular, Timesteps }; ReadingMode current_reading_mode_; boost::gregorian::date current_reading_date_; int current_epoch_; From 3aa067b3742c3f04df8a6d61aa4c93d6c99d3920 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 5 Jun 2012 14:06:00 +0200 Subject: [PATCH 04/23] Added WELOPEN, now schedules DATES and TSTEP properly. --- opm/core/eclipse/EclipseGridParser.cpp | 52 +++++++++++++++---- opm/core/eclipse/EclipseGridParser.hpp | 4 +- opm/core/eclipse/SpecialEclipseFields.hpp | 61 +++++++++++++++++++++++ 3 files changed, 107 insertions(+), 10 deletions(-) diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index 0780faaf..a91629c2 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -41,6 +41,7 @@ #include #include #include +#include #include #include #include @@ -97,6 +98,7 @@ namespace EclipseKeywords string("SGOF"), string("SWOF"), string("ROCK"), string("ROCKTAB"), string("WELSPECS"), string("COMPDAT"), string("WCONINJE"), string("WCONPROD"), string("WELTARG"), + string("WELOPEN"), string("EQUIL"), string("PVCDO"), string("TSTEP"), string("PLYVISC"), string("PLYROCK"), string("PLYADS"), string("PLYMAX"), string("TLMIXPAR"), string("WPOLYMER"), @@ -280,6 +282,8 @@ namespace { EclipseGridParser::EclipseGridParser() //--------------------------------------------------------------------------- : current_reading_mode_(Regular), + start_date_(boost::date_time::not_a_date_time), + current_time_days_(0.0), current_epoch_(0) { } @@ -290,6 +294,8 @@ EclipseGridParser::EclipseGridParser() EclipseGridParser::EclipseGridParser(const string& filename, bool convert_to_SI) //--------------------------------------------------------------------------- : current_reading_mode_(Regular), + start_date_(boost::date_time::not_a_date_time), + current_time_days_(0.0), current_epoch_(0) { // Store directory of filename @@ -384,22 +390,50 @@ void EclipseGridParser::readImpl(istream& is) break; } case Timestepping: { + SpecialMap& sm = special_field_by_epoch_[current_epoch_]; + if (start_date_.is_not_a_date()) { + // Set it to START date, or default if no START. + // This will only ever happen in the first epoch, + // upon first encountering a timestepping keyword. + SpecialMap::const_iterator it = sm.find("START"); + if (hasField("START")) { + start_date_ = getSTART().date; + } else { + start_date_ = boost::gregorian::date(1983, 1, 1); + } + } if (current_reading_mode_ == Regular) { current_reading_mode_ = Timesteps; } - SpecialMap& sm = special_field_by_epoch_[current_epoch_]; - // Append to current epoch's TSTEP, if it exists. + // Get current epoch's TSTEP, if it exists, create new if not. SpecialMap::iterator it = sm.find("TSTEP"); + TSTEP* tstep = 0; if (it != sm.end()) { - it->second->read(is); // This will append to the TSTEP object. + tstep = dynamic_cast(it->second.get()); } else { - // There is no existing TSTEP for this epoch, create it. - SpecialFieldPtr sb_ptr = createSpecialField(is, keyword); - if (sb_ptr) { - sm[keyword] = sb_ptr; - } else { - THROW("Could not create field " << keyword); + SpecialFieldPtr sb_ptr(new TSTEP()); + tstep = dynamic_cast(sb_ptr.get()); + sm["TSTEP"] = sb_ptr; + } + ASSERT(tstep != 0); + // Append new steps to current TSTEP object + if (keyword == "TSTEP") { + const int num_steps_old = tstep->tstep_.size(); + tstep->read(is); // This will append to the TSTEP object. + const double added_days + = std::accumulate(tstep->tstep_.begin() + num_steps_old, tstep->tstep_.end(), 0.0); + current_time_days_ += added_days; + } else if (keyword == "DATES") { + DATES dates; + dates.read(is); + for (std::size_t dix = 0; dix < dates.dates.size(); ++dix) { + boost::gregorian::date_duration since_start = dates.dates[dix] - start_date_; + double step = double(since_start.days()) - current_time_days_; + tstep->tstep_.push_back(step); + current_time_days_ = double(since_start.days()); } + } else { + THROW("Keyword " << keyword << " cannot be handled here."); } break; } diff --git a/opm/core/eclipse/EclipseGridParser.hpp b/opm/core/eclipse/EclipseGridParser.hpp index cf57906b..659b3484 100644 --- a/opm/core/eclipse/EclipseGridParser.hpp +++ b/opm/core/eclipse/EclipseGridParser.hpp @@ -150,6 +150,7 @@ public: SPECIAL_FIELD(WCONINJE); SPECIAL_FIELD(WCONPROD); SPECIAL_FIELD(WELTARG); + SPECIAL_FIELD(WELOPEN); SPECIAL_FIELD(EQUIL); SPECIAL_FIELD(PVCDO); SPECIAL_FIELD(TSTEP); @@ -205,7 +206,8 @@ private: // For SCHEDULE handling. enum ReadingMode { Regular, Timesteps }; ReadingMode current_reading_mode_; - boost::gregorian::date current_reading_date_; + boost::gregorian::date start_date_; + double current_time_days_; int current_epoch_; typedef std::map SpecialMap; std::vector special_field_by_epoch_; diff --git a/opm/core/eclipse/SpecialEclipseFields.hpp b/opm/core/eclipse/SpecialEclipseFields.hpp index d7e3548f..7e2fce9a 100644 --- a/opm/core/eclipse/SpecialEclipseFields.hpp +++ b/opm/core/eclipse/SpecialEclipseFields.hpp @@ -1658,6 +1658,67 @@ struct WELTARG : public SpecialBase }; + +struct WelopenLine +{ + std::string well_; // Well name or well name root + std::string openshutflag_; // What to do with the well. +}; + +/// Class for keyword WELOPEN +struct WELOPEN : public SpecialBase +{ + std::vector welopen; + + WELOPEN() + { + } + + virtual ~WELOPEN() + {} + + virtual std::string name() const {return std::string("WELOPEN");} + + virtual void read(std::istream& is) + { + while(is) { + std::string wellname = readString(is); + if (wellname[0] == '/') { + is >> ignoreLine; + break; + } + while (wellname.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + wellname = readString(is); + } + WelopenLine welopen_line; + welopen_line.well_ = wellname; + welopen_line.openshutflag_ = readString(is); + ignoreSlashLine(is); + welopen.push_back(welopen_line); + } + } + + virtual void write(std::ostream& os) const + { + os << name() << std::endl; + for (int i=0; i<(int)welopen.size(); ++i) { + os << welopen[i].well_ << " " + << welopen[i].openshutflag_ + << std::endl; + } + os << std::endl; + } + + virtual void convertToSI(const EclipseUnits& /*units*/) + { + } +}; + + + + /// Class holding a data line of keyword EQUIL struct EquilLine { From 831b2bb40bf558ce37271ef923846647bbb9a9b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 5 Jun 2012 15:20:13 +0200 Subject: [PATCH 05/23] Added skeleton doxygen_main.hpp for overview documentation. Also add generation of docs from opm/core, and disable docs from unfinished tutorial. --- Doxyfile | 4 ++-- opm/core/doxygen_main.hpp | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 37 insertions(+), 2 deletions(-) create mode 100644 opm/core/doxygen_main.hpp diff --git a/Doxyfile b/Doxyfile index 4b1c43ee..3f1fb256 100644 --- a/Doxyfile +++ b/Doxyfile @@ -610,8 +610,8 @@ WARN_LOGFILE = # directories like "/usr/src/myproject". Separate the files or directories # with spaces. -#INPUT = opm/core tutorials examples -INPUT = tutorials +INPUT = opm/core tutorials/tutorial1.cpp tutorials/tutorial2.cpp tutorials/tutorial3.cpp examples +#INPUT = tutorials # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is diff --git a/opm/core/doxygen_main.hpp b/opm/core/doxygen_main.hpp new file mode 100644 index 00000000..e33af59a --- /dev/null +++ b/opm/core/doxygen_main.hpp @@ -0,0 +1,35 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_DOXYGEN_MAIN_HEADER_INCLUDED +#define OPM_DOXYGEN_MAIN_HEADER_INCLUDED + + +/** \mainpage Documentation for the opm-core library. + +The following are the main library features: +

Grid handling

+

Well handling

+

Pressure solvers

+

Transport solvers

+

Various utilities

+ +*/ + +#endif // OPM_DOXYGEN_MAIN_HEADER_INCLUDED From 955b5ea06811c723b42b672acb7b74e3b639be26 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 5 Jun 2012 15:42:49 +0200 Subject: [PATCH 06/23] Reorganized, added opm/core/wells/ and opm/core/simulator/. --- Makefile.am | 32 +++++++++---------- examples/sim_wateroil.cpp | 8 ++--- examples/spu_2p.cpp | 6 ++-- examples/wells_example.cpp | 4 +-- opm/core/{utility => }/newwells.c | 0 opm/core/pressure/CompressibleTpfa.cpp | 4 +-- opm/core/{ => simulator}/BlackoilState.hpp | 0 opm/core/{ => simulator}/TwophaseState.hpp | 0 opm/core/{ => simulator}/WellState.hpp | 0 opm/core/{ => utility}/ColumnExtract.hpp | 0 .../{ => wells}/InjectionSpecification.cpp | 3 +- .../{ => wells}/InjectionSpecification.hpp | 0 .../{ => wells}/ProductionSpecification.cpp | 2 +- .../{ => wells}/ProductionSpecification.hpp | 0 opm/core/{ => wells}/WellCollection.cpp | 2 +- opm/core/{ => wells}/WellCollection.hpp | 2 +- opm/core/{ => wells}/WellsGroup.cpp | 2 +- opm/core/{ => wells}/WellsGroup.hpp | 4 +-- opm/core/{ => wells}/WellsManager.cpp | 4 +-- opm/core/{ => wells}/WellsManager.hpp | 4 +-- tests/test_column_extract.cpp | 2 +- tutorials/tutorial3.cpp | 2 +- tutorials/tutorial4.cpp | 4 +-- 23 files changed, 43 insertions(+), 42 deletions(-) rename opm/core/{utility => }/newwells.c (100%) rename opm/core/{ => simulator}/BlackoilState.hpp (100%) rename opm/core/{ => simulator}/TwophaseState.hpp (100%) rename opm/core/{ => simulator}/WellState.hpp (100%) rename opm/core/{ => utility}/ColumnExtract.hpp (100%) rename opm/core/{ => wells}/InjectionSpecification.cpp (88%) rename opm/core/{ => wells}/InjectionSpecification.hpp (100%) rename opm/core/{ => wells}/ProductionSpecification.cpp (88%) rename opm/core/{ => wells}/ProductionSpecification.hpp (100%) rename opm/core/{ => wells}/WellCollection.cpp (99%) rename opm/core/{ => wells}/WellCollection.hpp (99%) rename opm/core/{ => wells}/WellsGroup.cpp (99%) rename opm/core/{ => wells}/WellsGroup.hpp (99%) rename opm/core/{ => wells}/WellsManager.cpp (99%) rename opm/core/{ => wells}/WellsManager.hpp (98%) diff --git a/Makefile.am b/Makefile.am index a445e380..db4a48f5 100644 --- a/Makefile.am +++ b/Makefile.am @@ -51,6 +51,8 @@ opm/core/fluid/RockFromDeck.cpp \ opm/core/fluid/SaturationPropsBasic.cpp \ opm/core/fluid/SaturationPropsFromDeck.cpp \ opm/core/grid.c \ +opm/core/newwells.c \ +opm/core/GridManager.cpp \ opm/core/grid/cart_grid.c \ opm/core/grid/cornerpoint_grid.c \ opm/core/grid/cpgpreprocess/geometry.c \ @@ -72,14 +74,12 @@ opm/core/utility/miscUtilities.cpp \ opm/core/utility/miscUtilitiesBlackoil.cpp \ opm/core/utility/SimulatorTimer.cpp \ opm/core/utility/StopWatch.cpp \ -opm/core/utility/newwells.c \ opm/core/utility/writeVtkData.cpp \ -opm/core/GridManager.cpp \ -opm/core/WellsManager.cpp \ -opm/core/WellsGroup.cpp \ -opm/core/WellCollection.cpp \ -opm/core/InjectionSpecification.cpp \ -opm/core/ProductionSpecification.cpp \ +opm/core/wells/WellsManager.cpp \ +opm/core/wells/WellsGroup.cpp \ +opm/core/wells/WellCollection.cpp \ +opm/core/wells/InjectionSpecification.cpp \ +opm/core/wells/ProductionSpecification.cpp \ opm/core/linalg/sparse_sys.c \ opm/core/linalg/LinearSolverInterface.cpp \ opm/core/linalg/LinearSolverFactory.cpp \ @@ -158,6 +158,7 @@ opm/core/grid/cpgpreprocess/geometry.h \ opm/core/grid/cpgpreprocess/facetopology.h \ opm/core/grid/cpgpreprocess/grdecl.h \ opm/core/utility/Average.hpp \ +opm/core/utility/ColumnExtract.hpp \ opm/core/utility/ErrorMacros.hpp \ opm/core/utility/Factory.hpp \ opm/core/utility/MonotCubicInterpolator.hpp \ @@ -193,17 +194,16 @@ opm/core/linalg/LinearSolverFactory.hpp \ opm/core/grid.h \ opm/core/well.h \ opm/core/newwells.h \ -opm/core/WellsGroup.hpp \ -opm/core/WellCollection.hpp \ -opm/core/InjectionSpecification.hpp \ -opm/core/ProductionSpecification.hpp \ -opm/core/BlackoilState.hpp \ -opm/core/ColumnExtract.hpp \ +opm/core/wells/WellsManager.hpp \ +opm/core/wells/WellsGroup.hpp \ +opm/core/wells/WellCollection.hpp \ +opm/core/wells/InjectionSpecification.hpp \ +opm/core/wells/ProductionSpecification.hpp \ opm/core/GridAdapter.hpp \ opm/core/GridManager.hpp \ -opm/core/TwophaseState.hpp \ -opm/core/WellsManager.hpp \ -opm/core/WellState.hpp \ +opm/core/simulator/BlackoilState.hpp \ +opm/core/simulator/TwophaseState.hpp \ +opm/core/simulator/WellState.hpp \ opm/core/pressure/fsh.h \ opm/core/pressure/HybridPressureSolver.hpp \ opm/core/pressure/IncompTpfa.hpp \ diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp index 6c26c0fa..e2421bb7 100644 --- a/examples/sim_wateroil.cpp +++ b/examples/sim_wateroil.cpp @@ -27,7 +27,7 @@ #include #include #include -#include +#include #include #include #include @@ -44,9 +44,9 @@ #include -#include -#include -#include +#include +#include +#include #include #include diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 30d6f71e..490d9648 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -43,7 +43,7 @@ #include #include #include -#include +#include #include #include #include @@ -69,8 +69,8 @@ #include #include -#include -#include +#include +#include #include #include diff --git a/examples/wells_example.cpp b/examples/wells_example.cpp index 18c051aa..a8bdaae9 100644 --- a/examples/wells_example.cpp +++ b/examples/wells_example.cpp @@ -5,14 +5,14 @@ #include "opm/core/utility/initState.hpp" #include "opm/core/utility/SimulatorTimer.hpp" -#include +#include #include #include #include #include #include #include -#include +#include #include #include #include diff --git a/opm/core/utility/newwells.c b/opm/core/newwells.c similarity index 100% rename from opm/core/utility/newwells.c rename to opm/core/newwells.c diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 0ee8bf8b..83d3b804 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -28,8 +28,8 @@ #include #include #include -#include -#include +#include +#include #include #include diff --git a/opm/core/BlackoilState.hpp b/opm/core/simulator/BlackoilState.hpp similarity index 100% rename from opm/core/BlackoilState.hpp rename to opm/core/simulator/BlackoilState.hpp diff --git a/opm/core/TwophaseState.hpp b/opm/core/simulator/TwophaseState.hpp similarity index 100% rename from opm/core/TwophaseState.hpp rename to opm/core/simulator/TwophaseState.hpp diff --git a/opm/core/WellState.hpp b/opm/core/simulator/WellState.hpp similarity index 100% rename from opm/core/WellState.hpp rename to opm/core/simulator/WellState.hpp diff --git a/opm/core/ColumnExtract.hpp b/opm/core/utility/ColumnExtract.hpp similarity index 100% rename from opm/core/ColumnExtract.hpp rename to opm/core/utility/ColumnExtract.hpp diff --git a/opm/core/InjectionSpecification.cpp b/opm/core/wells/InjectionSpecification.cpp similarity index 88% rename from opm/core/InjectionSpecification.cpp rename to opm/core/wells/InjectionSpecification.cpp index dca7fc76..accae78c 100644 --- a/opm/core/InjectionSpecification.cpp +++ b/opm/core/wells/InjectionSpecification.cpp @@ -1,4 +1,5 @@ -#include +#include + namespace Opm { diff --git a/opm/core/InjectionSpecification.hpp b/opm/core/wells/InjectionSpecification.hpp similarity index 100% rename from opm/core/InjectionSpecification.hpp rename to opm/core/wells/InjectionSpecification.hpp diff --git a/opm/core/ProductionSpecification.cpp b/opm/core/wells/ProductionSpecification.cpp similarity index 88% rename from opm/core/ProductionSpecification.cpp rename to opm/core/wells/ProductionSpecification.cpp index eb0f7d83..4d4bdf4d 100644 --- a/opm/core/ProductionSpecification.cpp +++ b/opm/core/wells/ProductionSpecification.cpp @@ -1,4 +1,4 @@ -#include +#include namespace Opm { diff --git a/opm/core/ProductionSpecification.hpp b/opm/core/wells/ProductionSpecification.hpp similarity index 100% rename from opm/core/ProductionSpecification.hpp rename to opm/core/wells/ProductionSpecification.hpp diff --git a/opm/core/WellCollection.cpp b/opm/core/wells/WellCollection.cpp similarity index 99% rename from opm/core/WellCollection.cpp rename to opm/core/wells/WellCollection.cpp index 76b311e7..727a7242 100644 --- a/opm/core/WellCollection.cpp +++ b/opm/core/wells/WellCollection.cpp @@ -18,7 +18,7 @@ You should have received a copy of the GNU General Public License along with OpenRS. If not, see . */ -#include "WellCollection.hpp" +#include namespace Opm { diff --git a/opm/core/WellCollection.hpp b/opm/core/wells/WellCollection.hpp similarity index 99% rename from opm/core/WellCollection.hpp rename to opm/core/wells/WellCollection.hpp index 62abd4ee..1957576a 100644 --- a/opm/core/WellCollection.hpp +++ b/opm/core/wells/WellCollection.hpp @@ -23,7 +23,7 @@ along with OpenRS. If not, see . #define OPM_WELLCOLLECTION_HPP #include -#include +#include #include #include diff --git a/opm/core/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp similarity index 99% rename from opm/core/WellsGroup.cpp rename to opm/core/wells/WellsGroup.cpp index 903cc9b2..3514bb5d 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -17,7 +17,7 @@ along with OPM. If not, see . */ -#include +#include #include #include #include diff --git a/opm/core/WellsGroup.hpp b/opm/core/wells/WellsGroup.hpp similarity index 99% rename from opm/core/WellsGroup.hpp rename to opm/core/wells/WellsGroup.hpp index 3c14950d..f1b1afbb 100644 --- a/opm/core/WellsGroup.hpp +++ b/opm/core/wells/WellsGroup.hpp @@ -20,8 +20,8 @@ #ifndef OPM_WELLSGROUP_HPP #define OPM_WELLSGROUP_HPP -#include -#include +#include +#include #include #include #include diff --git a/opm/core/WellsManager.cpp b/opm/core/wells/WellsManager.cpp similarity index 99% rename from opm/core/WellsManager.cpp rename to opm/core/wells/WellsManager.cpp index 8aa5d644..d305ee9f 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -18,13 +18,13 @@ */ -#include +#include #include #include #include #include #include -#include +#include #include #include diff --git a/opm/core/WellsManager.hpp b/opm/core/wells/WellsManager.hpp similarity index 98% rename from opm/core/WellsManager.hpp rename to opm/core/wells/WellsManager.hpp index 3fbdbd1a..6df65ec4 100644 --- a/opm/core/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -21,8 +21,8 @@ #define OPM_WELLSMANAGER_HEADER_INCLUDED -#include -#include +#include +#include struct Wells; struct UnstructuredGrid; diff --git a/tests/test_column_extract.cpp b/tests/test_column_extract.cpp index 877f21a0..46b1d538 100644 --- a/tests/test_column_extract.cpp +++ b/tests/test_column_extract.cpp @@ -6,7 +6,7 @@ #define BOOST_TEST_MODULE ColumnExtractTest #include -#include +#include #include #include diff --git a/tutorials/tutorial3.cpp b/tutorials/tutorial3.cpp index 6420c199..fc15ab30 100644 --- a/tutorials/tutorial3.cpp +++ b/tutorials/tutorial3.cpp @@ -37,7 +37,7 @@ #include -#include +#include #include #include diff --git a/tutorials/tutorial4.cpp b/tutorials/tutorial4.cpp index 0897b5b4..62776eaf 100644 --- a/tutorials/tutorial4.cpp +++ b/tutorials/tutorial4.cpp @@ -44,12 +44,12 @@ #include #include -#include +#include #include #include #include -#include +#include /// \page tutorial4 Well controls /// From 40ef0e9573d3ba4b0893562414370c2ba41c9fa3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 5 Jun 2012 18:40:36 +0200 Subject: [PATCH 07/23] Revamp Fortran support in preparation of enabling Notay's AGMG solver. Specifically: - Enable C++, Fortran 77 and Fortran (95+) through Libtool's LT_LANG if available and through the traditional AC_PROG_* macros if not. This configuration is compatible with the versions of Libtool easily available for testing. For whatever reason--possibly a programming error in Libtool proper--invoking the AC_PROG_F* macros either directly or through AC_REQUIRE following an LT_INIT invocation leads to various ``expanded before required'' warnings. Searching the Autotools mailing lists does suggest that the interaction of C++, Fortran and Libtool is traditionally somewhat unstable but has improved in very recent editions of Autoconf and Libtool. - Re-factor the LAPACK support out to a custom macro, OPM_LAPACK, and invoke it from configure.ac. --- configure.ac | 25 +++++++++++++------------ m4/opm_lapack.m4 | 4 ++++ 2 files changed, 17 insertions(+), 12 deletions(-) create mode 100644 m4/opm_lapack.m4 diff --git a/configure.ac b/configure.ac index 5b85b794..ff1d6d95 100644 --- a/configure.ac +++ b/configure.ac @@ -8,32 +8,33 @@ AM_INIT_AUTOMAKE([-Wall -Werror foreign subdir-objects]) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) -m4_ifdef([LT_INIT], [LT_INIT]) - AC_CONFIG_MACRO_DIR([m4]) AC_CONFIG_SRCDIR([opm/core/grid.h]) AC_CONFIG_HEADERS([config.h]) # Checks for programs. -AC_PROG_CXX AC_PROG_CC AM_PROG_CC_C_O -m4_ifdef([LT_INIT], [], [AC_PROG_LIBTOOL]) -AC_PROG_RANLIB -# AX_LAPACK requires a working F77 compiler or, rather, its runtime -# support libraries. -AC_PROG_F77 - -# F77 name mangling -AC_F77_WRAPPERS +m4_ifdef([LT_INIT], + [LT_INIT[]dnl + LT_LANG([C++])dnl + LT_LANG([Fortran 77])dnl + LT_LANG([Fortran])dnl + ],dnl + [AC_PROG_LIBTOOL[]dnl + AC_PROG_CXX[]dnl + AC_PROG_F77[]dnl + AC_PROG_FC[]dnl + ])[]dnl # Checks for libraries. # Bring in numerics support (standard library component) AC_SEARCH_LIBS([sqrt], [m]) -AX_LAPACK +OPM_LAPACK + AX_BOOST_BASE([1.37]) AX_BOOST_SYSTEM AX_BOOST_DATE_TIME diff --git a/m4/opm_lapack.m4 b/m4/opm_lapack.m4 new file mode 100644 index 00000000..b5fa5ca5 --- /dev/null +++ b/m4/opm_lapack.m4 @@ -0,0 +1,4 @@ +AC_DEFUN([OPM_LAPACK], +[AC_REQUIRE([AC_F77_WRAPPERS])dnl + AC_REQUIRE([AX_LAPACK])dnl +])[]dnl From c724bb6b3a845937fd70738b1cb6c0b94f0172cb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 5 Jun 2012 19:44:18 +0200 Subject: [PATCH 08/23] Add preliminary build glue to include AGMG in the OPM-Core library. The approach will likely be changed due to the presence of "subdir-objects" in the AM_INIT_AUTOMAKE call. --- Makefile.am | 10 ++++++++++ configure.ac | 1 + m4/agmg.m4 | 29 +++++++++++++++++++++++++++++ 3 files changed, 40 insertions(+) create mode 100644 m4/agmg.m4 diff --git a/Makefile.am b/Makefile.am index db4a48f5..0b47f16c 100644 --- a/Makefile.am +++ b/Makefile.am @@ -271,3 +271,13 @@ opm/core/linalg/LinearSolverIstl.cpp nobase_include_HEADERS += \ opm/core/linalg/LinearSolverIstl.hpp endif + + +if BUILD_AGMG +libopmcore_la_SOURCES += \ +$(AGMG_SRCDIR)/dagmg.f90 \ +$(AGMG_SRCDIR)/dagmg_mumps.f90 + +libopmcore_la_LDFLAGS += \ +$(FCLIBS) +endif diff --git a/configure.ac b/configure.ac index ff1d6d95..7f82ebcc 100644 --- a/configure.ac +++ b/configure.ac @@ -42,6 +42,7 @@ AX_BOOST_FILESYSTEM AX_BOOST_UNIT_TEST_FRAMEWORK AX_DUNE_ISTL +OPM_AGMG # Checks for header files. AC_CHECK_HEADERS([float.h limits.h stddef.h stdlib.h string.h]) diff --git a/m4/agmg.m4 b/m4/agmg.m4 new file mode 100644 index 00000000..d8d446e6 --- /dev/null +++ b/m4/agmg.m4 @@ -0,0 +1,29 @@ +AC_DEFUN([OPM_F90_COMPILATION_SYSTEM], +[AC_REQUIRE([AC_PROG_FC_C_O])dnl + AC_REQUIRE([AC_FC_WRAPPERS])dnl + AC_REQUIRE([AC_FC_LIBRARY_LDFLAGS])dnl +])[]dnl + +AC_DEFUN([OPM_AGMG],dnl +[AC_ARG_WITH([agmg],dnl + [AS_HELP_STRING([--with-agmg=SRCDIR],dnl + [Include DOUBLE PRECISION version Notay's of AGMG Algebraic + Multigrid solver from specified source location SRCDIR. Note: this + option requires a complete, working Fortran 90 environment.])], + [],dnl + [with_agmg=no])[]dnl + + AS_IF([test x"$with_agmg" != x"no" -a -d "$with_agmg"],dnl + [AS_IF([test -f "$with_agmg/dagmg.f90"],dnl + [AC_SUBST([AGMG_SRCDIR], [$with_agmg])[]dnl + AC_DEFINE([HAVE_AGMG], [1],dnl + [Define to `1' if Notay's AGMG solver is included])[]dnl + build_agmg="yes"],dnl + [build_agmg="no"])],dnl + [build_agmg="no"])[]dnl + + AS_IF([test x"$build_agmg" = x"yes"],dnl + [AC_REQUIRE([OPM_F90_COMPILATION_SYSTEM])], [:])[]dnl + + AM_CONDITIONAL([BUILD_AGMG], [test x"$build_agmg" = x"yes"]) +])[]dnl From 496be60220f19b4115e2e15069749f3bbf06bb75 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 5 Jun 2012 19:56:17 +0200 Subject: [PATCH 09/23] Don't use AC_REQUIRE within a conditional. Older versions of Autoconf are not able to place the tests at appropriate locations. --- m4/agmg.m4 | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/m4/agmg.m4 b/m4/agmg.m4 index d8d446e6..9e09157b 100644 --- a/m4/agmg.m4 +++ b/m4/agmg.m4 @@ -1,9 +1,3 @@ -AC_DEFUN([OPM_F90_COMPILATION_SYSTEM], -[AC_REQUIRE([AC_PROG_FC_C_O])dnl - AC_REQUIRE([AC_FC_WRAPPERS])dnl - AC_REQUIRE([AC_FC_LIBRARY_LDFLAGS])dnl -])[]dnl - AC_DEFUN([OPM_AGMG],dnl [AC_ARG_WITH([agmg],dnl [AS_HELP_STRING([--with-agmg=SRCDIR],dnl @@ -23,7 +17,10 @@ AC_DEFUN([OPM_AGMG],dnl [build_agmg="no"])[]dnl AS_IF([test x"$build_agmg" = x"yes"],dnl - [AC_REQUIRE([OPM_F90_COMPILATION_SYSTEM])], [:])[]dnl + [AC_PROG_FC_C_O[]dnl + AC_FC_WRAPPERS[]dnl + AC_FC_LIBRARY_LDFLAGS[]dnl + ], [:])[]dnl AM_CONDITIONAL([BUILD_AGMG], [test x"$build_agmg" = x"yes"]) ])[]dnl From ca7b4a5da3cc9db09aadbe3c14b1f174f8fbbd77 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 00:26:10 +0200 Subject: [PATCH 10/23] Add simple gateway to Notay's AGMG solver package. Hook up to build, but don't expose through the LinearSolverFactory. Only build tested. Won't be built unless the user explicitly enables AGMG support. --- Makefile.am | 6 +- opm/core/linalg/LinearSolverAGMG.cpp | 122 +++++++++++++++++++++++++++ opm/core/linalg/LinearSolverAGMG.hpp | 96 +++++++++++++++++++++ 3 files changed, 223 insertions(+), 1 deletion(-) create mode 100644 opm/core/linalg/LinearSolverAGMG.cpp create mode 100644 opm/core/linalg/LinearSolverAGMG.hpp diff --git a/Makefile.am b/Makefile.am index 0b47f16c..c64ada2a 100644 --- a/Makefile.am +++ b/Makefile.am @@ -276,7 +276,11 @@ endif if BUILD_AGMG libopmcore_la_SOURCES += \ $(AGMG_SRCDIR)/dagmg.f90 \ -$(AGMG_SRCDIR)/dagmg_mumps.f90 +$(AGMG_SRCDIR)/dagmg_mumps.f90 \ +opm/core/linalg/LinearSolverAGMG.cpp + +nobase_include_HEADERS += \ +opm/core/linalg/LinearSolverAGMG.hpp libopmcore_la_LDFLAGS += \ $(FCLIBS) diff --git a/opm/core/linalg/LinearSolverAGMG.cpp b/opm/core/linalg/LinearSolverAGMG.cpp new file mode 100644 index 00000000..d7b47089 --- /dev/null +++ b/opm/core/linalg/LinearSolverAGMG.cpp @@ -0,0 +1,122 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#if HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include +#include +#include +#include +#include + +#if HAVE_AGMG +// Manual prototype of main gateway routine to DOUBLE PRECISION, +// serial version of the AGMG software package. +// +// Note that both the matrix entries and column indices are writable. +// The solver may permute the matrix entries within each row during +// the setup phase. +#define DAGMG_ F77_FUNC(dagmg, DAGMG) + +extern "C" +void +DAGMG_(const int* N , // System size + double* sa , // Non-zero elements + int* ja , // Column indices + const int* ia , // Row pointers + const double* f , // Right-hand side + double* x , // Solution + int* ijob , // Main control parameter + int* iprint, // Message output unit + int* nrest , // Number of GCR restarts + int* iter , // Maximum (and actual) number of iterations + const double* tol ); // Residual reduction tolerance + +#endif // HAVE_AGMG + +namespace Opm +{ + LinearSolverAGMG::LinearSolverAGMG(const int max_it, + const double rtol , + const bool is_spd) + : max_it_(max_it), + rtol_ (rtol) , + is_spd_(is_spd) + { +#if !HAVE_AGMG + THROW("AGMG support is not enabled in this library"); +#endif // HAVE_AGMG + } + + LinearSolverInterface::LinearSolverReport + LinearSolverAGMG::solve(const int size , + const int nonzeros, + const int* ia , + const int* ja , + const double* sa , + const double* rhs , + double* solution) const + { + const std::vector::size_type nnz = ia[size]; + + assert (nnz == std::vector::size_type(nonzeros)); + +#if defined(NDEBUG) + // Suppress warning about unused parameter. + static_cast(nonzeros); +#endif + + std::vector a(sa, sa + nnz); + + // Account for 1-based indexing. + std::vector i(ia, ia + std::vector::size_type(size + 1)); + std::transform(i.begin(), i.end(), i.begin(), + std::bind2nd(std::plus(), 1)); + + std::vector j(ja, ja + nnz); + std::transform(j.begin(), j.end(), j.begin(), + std::bind2nd(std::plus(), 1)); + + LinearSolverInterface::LinearSolverReport rpt = {}; + rpt.iterations = max_it_; + + int ijob = 0; // Setup + solution + cleanup, x0==0. + + int nrest; + if (is_spd_) { + nrest = 1; // Use CG algorithm + } + else { + nrest = 10; // Suggested default number of GCR restarts. + } + + int iprint = 0; // Suppress most output + DAGMG_(& size, & a[0], & j[0], & i[0], rhs, solution, + & ijob, & iprint, & nrest, & rpt.iterations, & rtol_); + + rpt.converged = rpt.iterations <= max_it_; + return rpt; + } +} // namespace Opm diff --git a/opm/core/linalg/LinearSolverAGMG.hpp b/opm/core/linalg/LinearSolverAGMG.hpp new file mode 100644 index 00000000..d53debba --- /dev/null +++ b/opm/core/linalg/LinearSolverAGMG.hpp @@ -0,0 +1,96 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_LINEARSOLVERAGMG_HEADER_INCLUDED +#define OPM_LINEARSOLVERAGMG_HEADER_INCLUDED + +/** + * \file + * Gateway to Notay's AGMG package implementing an aggregation-based + * algebraic multigrid method. + * + * References: + * * Y. Notay, + * An aggregation-based algebraic multigrid method, + * Electronic Transactions on Numerical Analysis, vol 37, + * pp. 123-146, 2010. + * + * * A. Napov and Y. Notay, + * An algebraic multigrid method with guaranteed convergence rate, + * Report GANMN 10-03, Université Libre de Bruxelles, Brussels, + * Belgium, 2010 (Revised 2011). + * + * * Y. Notay, + * Aggregation-based algebraic multigrid for convection-diffusion + * equations, Report GANMN 11-01, Université Libre de Bruxelles, + * Brussels, Belgium, 2011. + */ + +#include + +namespace Opm +{ + /// Concrete class encapsulating Notay's AGMG package + class LinearSolverAGMG : public LinearSolverInterface + { + public: + /** + * Constructor. + * \param[in] max_it Maximum number of solver iterations. + * \param[in] rtol Residual reduction tolerance. + * \param[in] is_spd Whether or not the matrix is SPD. SPD + * systems are solved using CG while other + * systems are solved using GCR. + */ + LinearSolverAGMG(const int max_it = 100 , + const double rtol = 1.0e-6, + const bool is_spd = false); + + /// Destructor. + virtual ~LinearSolverAGMG(); + + using LinearSolverInterface::solve; + + /// Solve a linear system, with a matrix given in compressed + /// sparse row format. + /// \param[in] size Number of rows (and colums). + /// \param[in] nonzeros Number of (structural) non-zeros. + /// \param[in] ia Row pointers. + /// \param[in] ja Column indices. + /// \param[in] sa (structurally) non-zero elements. + /// \param[in] rhs System right-hand side. + /// \param[inout] solution System solution. + /// \return Solver meta-data concerning most recent system solve. + virtual LinearSolverInterface::LinearSolverReport + solve(const int size, const int nonzeros, + const int* ia, const int* ja, const double* sa, + const double* rhs, double* solution) const; + + private: + int max_it_; + double rtol_ ; + bool is_spd_; + }; + + +} // namespace Opm + + + +#endif // OPM_LINEARSOLVERISTL_HEADER_INCLUDED From 61a5f3d1e54f47fbe245f3db37782ebbdadb2bdb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 00:27:10 +0200 Subject: [PATCH 11/23] Fix misprint left over from early development. --- opm/core/linalg/LinearSolverAGMG.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/linalg/LinearSolverAGMG.cpp b/opm/core/linalg/LinearSolverAGMG.cpp index d7b47089..c6716f5c 100644 --- a/opm/core/linalg/LinearSolverAGMG.cpp +++ b/opm/core/linalg/LinearSolverAGMG.cpp @@ -38,7 +38,7 @@ // Note that both the matrix entries and column indices are writable. // The solver may permute the matrix entries within each row during // the setup phase. -#define DAGMG_ F77_FUNC(dagmg, DAGMG) +#define DAGMG_ FC_FUNC(dagmg, DAGMG) extern "C" void From 79c437ab405a44e29e00f838047b04e9a49f97b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 00:30:12 +0200 Subject: [PATCH 12/23] Fix another misprint from early implementation. --- opm/core/linalg/LinearSolverAGMG.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/linalg/LinearSolverAGMG.hpp b/opm/core/linalg/LinearSolverAGMG.hpp index d53debba..6056dc90 100644 --- a/opm/core/linalg/LinearSolverAGMG.hpp +++ b/opm/core/linalg/LinearSolverAGMG.hpp @@ -93,4 +93,4 @@ namespace Opm -#endif // OPM_LINEARSOLVERISTL_HEADER_INCLUDED +#endif // OPM_LINEARSOLVERAGMG_HEADER_INCLUDED From 66cbaf01b3ce8b4c882ac3603a0432bfe7f1cd2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 6 Jun 2012 09:12:15 +0200 Subject: [PATCH 13/23] Work in progress on simulator class for two-phase simulation. --- Makefile.am | 8 +- opm/core/simulator/SimulatorTwophase.cpp | 112 +++++++++++++++++++++++ opm/core/simulator/SimulatorTwophase.hpp | 62 +++++++++++++ 3 files changed, 179 insertions(+), 3 deletions(-) create mode 100644 opm/core/simulator/SimulatorTwophase.cpp create mode 100644 opm/core/simulator/SimulatorTwophase.hpp diff --git a/Makefile.am b/Makefile.am index db4a48f5..259f03bd 100644 --- a/Makefile.am +++ b/Makefile.am @@ -109,6 +109,7 @@ opm/core/pressure/ifsh.c \ opm/core/pressure/mimetic/mimetic.c \ opm/core/pressure/mimetic/hybsys_global.c \ opm/core/pressure/mimetic/hybsys.c \ +opm/core/simulator/SimulatorTwophase.cpp \ opm/core/transport/spu_explicit.c \ opm/core/transport/spu_implicit.c \ opm/core/transport/transport_source.c \ @@ -201,9 +202,6 @@ opm/core/wells/InjectionSpecification.hpp \ opm/core/wells/ProductionSpecification.hpp \ opm/core/GridAdapter.hpp \ opm/core/GridManager.hpp \ -opm/core/simulator/BlackoilState.hpp \ -opm/core/simulator/TwophaseState.hpp \ -opm/core/simulator/WellState.hpp \ opm/core/pressure/fsh.h \ opm/core/pressure/HybridPressureSolver.hpp \ opm/core/pressure/IncompTpfa.hpp \ @@ -230,6 +228,10 @@ opm/core/pressure/fsh_common_impl.h \ opm/core/pressure/mimetic/hybsys_global.h \ opm/core/pressure/mimetic/hybsys.h \ opm/core/pressure/mimetic/mimetic.h \ +opm/core/simulator/BlackoilState.hpp \ +opm/core/simulator/SimulatorTwophase.hpp \ +opm/core/simulator/TwophaseState.hpp \ +opm/core/simulator/WellState.hpp \ opm/core/transport/ImplicitTransport.hpp \ opm/core/transport/CSRMatrixUmfpackSolver.hpp \ opm/core/transport/CSRMatrixBlockAssembler.hpp \ diff --git a/opm/core/simulator/SimulatorTwophase.cpp b/opm/core/simulator/SimulatorTwophase.cpp new file mode 100644 index 00000000..a83b57b7 --- /dev/null +++ b/opm/core/simulator/SimulatorTwophase.cpp @@ -0,0 +1,112 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include + +namespace Opm +{ + + class SimulatorTwophase::Impl + { + public: + Impl() {} + void init(const parameter::ParameterGroup& param); + void run(const SimulatorTimer& timer, const Wells& wells, TwophaseState& state, WellState& well_state); + private: +#if 0 + bool output_; + std::string output_dir_; + int output_interval; + int num_transport_substeps_; + double gravity_[3]; + bool use_gravity_; + bool use_segregation_split = false; // + int nl_pressure_maxiter = 0; + double nl_pressure_tolerance = 0.0; + const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); + const int nl_maxiter = param.getDefault("nl_maxiter", 30); + + // boost::scoped_ptr grid_; + // boost::scoped_ptr props_; + // boost::scoped_ptr wells_; + // boost::scoped_ptr rock_comp_; + // Opm::SimulatorTimer simtimer_; + // Opm::TwophaseState state_; + // std::vector src(num_cells, 0.0); + // Opm::FlowBCManager bcs; + // Opm::LinearSolverFactory linsolver(param); + + + + + std::vector totmob; + std::vector omega; // Will remain empty if no gravity. + std::vector rc; // Will remain empty if no rock compressibility. + std::vector porevol; + double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + std::vector transport_src; + + Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver); + + Opm::TransportModelTwophase tsolver(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); + + typedef std::pair, std::vector > > ColMap; + ColMap columns; + + std::vector allcells(num_cells); +#endif + }; + + SimulatorTwophase::SimulatorTwophase() + { + pimpl_.reset(new Impl); + } + + + + void SimulatorTwophase::init(const parameter::ParameterGroup& param) + { + pimpl_->init(param); + } + + + + void SimulatorTwophase::run(const SimulatorTimer& timer, const Wells& wells, TwophaseState& state, WellState& well_state) + { + pimpl_->run(timer, wells, state, well_state); + } + + + + void SimulatorTwophase::Impl::init(const parameter::ParameterGroup& /*param*/) + { + } + + + + void SimulatorTwophase::Impl::run(const SimulatorTimer& timer, const Wells& wells, TwophaseState& state, WellState& well_state) + { + (void) timer; + (void) wells; + (void) state; + (void) well_state; + } + +} // namespace Opm diff --git a/opm/core/simulator/SimulatorTwophase.hpp b/opm/core/simulator/SimulatorTwophase.hpp new file mode 100644 index 00000000..636d4f94 --- /dev/null +++ b/opm/core/simulator/SimulatorTwophase.hpp @@ -0,0 +1,62 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_SIMULATORTWOPHASE_HEADER_INCLUDED +#define OPM_SIMULATORTWOPHASE_HEADER_INCLUDED + +#include + +struct UnstructuredGrid; +struct Wells; + +namespace Opm +{ + namespace parameter { class ParameterGroup; } + class SimulatorTimer; + class TwophaseState; + class WellState; + + /// Class collecting all necessary components for a two-phase simulation. + class SimulatorTwophase + { + public: + /// Default constructor. + SimulatorTwophase(); + + /// Initialise from parameters. + void init(const parameter::ParameterGroup& param); + + /// Run the simulation. + /// \param[in] timer governs the requested reporting timesteps + /// \param[in] wells data structure for wells + /// \param[out] state state of reservoir: pressure, fluxes + /// \param[out] well_state state of wells: bhp, perforation rates + void run(const SimulatorTimer& timer, + const Wells& wells, + TwophaseState& state, + WellState& well_state); + + private: + class Impl; + boost::scoped_ptr pimpl_; + }; + +} // namespace Opm + +#endif // OPM_SIMULATORTWOPHASE_HEADER_INCLUDED From 0a8eafc307a41466903a7faa65e203e20d95b9a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 09:46:17 +0200 Subject: [PATCH 14/23] Implement dtor to avoid linker errors. --- opm/core/linalg/LinearSolverAGMG.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/opm/core/linalg/LinearSolverAGMG.cpp b/opm/core/linalg/LinearSolverAGMG.cpp index c6716f5c..e028e2c5 100644 --- a/opm/core/linalg/LinearSolverAGMG.cpp +++ b/opm/core/linalg/LinearSolverAGMG.cpp @@ -70,6 +70,8 @@ namespace Opm #endif // HAVE_AGMG } + LinearSolverAGMG::~LinearSolverAGMG() {} + LinearSolverInterface::LinearSolverReport LinearSolverAGMG::solve(const int size , const int nonzeros, From 17ebf48f34ca1bea7af98c56a5a072ee3cc4ece3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 11:49:07 +0200 Subject: [PATCH 15/23] Add simple demonstration/regression test of AGMG solver. --- tests/Makefile.am | 6 +++ tests/test_agmg.cpp | 89 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+) create mode 100644 tests/test_agmg.cpp diff --git a/tests/Makefile.am b/tests/Makefile.am index a279e72d..064d385c 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -47,3 +47,9 @@ noinst_PROGRAMS += test_cfs_tpfa noinst_PROGRAMS += test_jacsys test_jacsys_SOURCES = test_jacsys.cpp endif + + +if BUILD_AGMG +noinst_PROGRAMS += test_agmg +test_agmg_SOURCES = test_agmg.cpp +endif diff --git a/tests/test_agmg.cpp b/tests/test_agmg.cpp new file mode 100644 index 00000000..b520468c --- /dev/null +++ b/tests/test_agmg.cpp @@ -0,0 +1,89 @@ +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace { + boost::shared_ptr + build_laplace_1d(const std::size_t m) + { + const std::size_t nnz = 3*(m - 2) + 2*2; + + boost::shared_ptr + A(csrmatrix_new_known_nnz(m, nnz), csrmatrix_delete); + + A->ia[ 0 ] = 0; + + // First row + A->ia[ 0 + 1 ] = A->ia[ 0 ]; + A->ja[ A->ia[0 + 1] ] = 0 + 0; + A->sa[ A->ia[0 + 1] ++ ] = 2.0; + A->ja[ A->ia[0 + 1] ] = 0 + 1; + A->sa[ A->ia[0 + 1] ++ ] = - 1.0; + + // General rows + for (std::size_t i = 1; i < m - 1; ++i) { + A->ia[i + 1] = A->ia[i]; + + A->ja[ A->ia[i + 1] ] = int(i) - 1; + A->sa[ A->ia[i + 1] ++ ] = - 1.0; + + A->ja[ A->ia[i + 1] ] = int(i) ; + A->sa[ A->ia[i + 1] ++ ] = 2.0; + + A->ja[ A->ia[i + 1] ] = int(i) + 1; + A->sa[ A->ia[i + 1] ++ ] = - 1.0; + } + + // Last row + A->ia[ (m - 1) + 1 ] = A->ia[ m - 1 ]; + + A->ja[ A->ia[ (m - 1) + 1 ] ] = int(m - 1) - 1; + A->sa[ A->ia[ (m - 1) + 1 ] ++ ] = - 1.0; + + A->ja[ A->ia[ (m - 1) + 1 ] ] = int(m - 1) ; + A->sa[ A->ia[ (m - 1) + 1 ] ++ ] = 2.0; + + return A; + } +} + + +int main() +{ + const std::size_t m = 10; + + boost::shared_ptr A = build_laplace_1d(m); + + // Form right-hand side [1, 0, 0, ...., 0, 1] + std::vector b(m, 0.0); + b[0] = 1.0; b.back() = 1.0; + + // Allocate solution vector + std::vector x(m); + + // Create solver for SPD system. + Opm::LinearSolverAGMG linsolve(100, 1e-9, true); + + Opm::LinearSolverInterface::LinearSolverReport + rpt = linsolve.solve(A.get(), & b[0], & x[0]); + + double e = 0.0; + for (std::size_t i = 0; i < m; ++i) { + const double d = x[i] - 1.0; + e += d * d; + } + + std::cerr << "|| e ||_2 = " + << std::scientific + << std::setprecision(5) + << std::sqrt(e) / double(m) << '\n'; + + return 0; +} From a6ce0197369e5c546e077fa0b61a900864f3545f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 12:38:27 +0200 Subject: [PATCH 16/23] Support matrices of two or more rows. --- tests/test_agmg.cpp | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/tests/test_agmg.cpp b/tests/test_agmg.cpp index b520468c..f04b3230 100644 --- a/tests/test_agmg.cpp +++ b/tests/test_agmg.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -10,10 +11,26 @@ #include namespace { + std::size_t + compute_nnz(const std::size_t m) + { + assert (m > 0); + + std::size_t nnz = m; // A(i,i) + + nnz += (m > 1) ? 2 : 0; // A(0,1), A(m-1,m-2) + nnz += (m > 2) ? 2 * (m - 2) : 0; // A(i,i-1), A(i,i+1) + + return nnz; + } + + boost::shared_ptr build_laplace_1d(const std::size_t m) { - const std::size_t nnz = 3*(m - 2) + 2*2; + assert (m >= 2); + + const std::size_t nnz = compute_nnz(m); boost::shared_ptr A(csrmatrix_new_known_nnz(m, nnz), csrmatrix_delete); From 47897c644cca1e42151c111401ccfee393618fa4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 12:39:29 +0200 Subject: [PATCH 17/23] Assert copyright. --- tests/test_agmg.cpp | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/tests/test_agmg.cpp b/tests/test_agmg.cpp index f04b3230..51c4a9cf 100644 --- a/tests/test_agmg.cpp +++ b/tests/test_agmg.cpp @@ -1,3 +1,23 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + Copyright 2012 Statoil ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + #include #include #include From 5702695879d3a47af9e2bcc961f81053c0c7c936 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 6 Jun 2012 13:42:25 +0200 Subject: [PATCH 18/23] Added some documentation. --- opm/core/simulator/WellState.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 56ede664..59d9f060 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -26,9 +26,11 @@ namespace Opm { + /// The state of a set of wells. class WellState { public: + /// Allocate and initialize if wells is non-null. template void init(const Wells* wells, const State& state) { @@ -44,9 +46,11 @@ namespace Opm } } + /// One bhp pressure per well. std::vector& bhp() { return bhp_; } const std::vector& bhp() const { return bhp_; } + /// One rate per well connection. std::vector& perfRates() { return perfrates_; } const std::vector& perfRates() const { return perfrates_; } From 3d054306278f71e20f5029f3bb9be428591ec0c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 6 Jun 2012 13:54:53 +0200 Subject: [PATCH 19/23] First version of SimulatorTwophase class done. Added test sim using it. --- examples/Makefile.am | 9 +- examples/sim_2p_incomp_reorder.cpp | 220 +++++++++++ opm/core/simulator/SimulatorTwophase.cpp | 471 ++++++++++++++++++++--- opm/core/simulator/SimulatorTwophase.hpp | 59 ++- 4 files changed, 684 insertions(+), 75 deletions(-) create mode 100644 examples/sim_2p_incomp_reorder.cpp diff --git a/examples/Makefile.am b/examples/Makefile.am index a60d05ea..7a178da6 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -8,7 +8,8 @@ noinst_PROGRAMS = \ refine_wells \ scaneclipsedeck \ sim_wateroil \ -wells_example +wells_example \ +sim_2p_incomp_reorder refine_wells_SOURCES = refine_wells.cpp @@ -18,6 +19,12 @@ $(LDADD) $(LIBS) \ $(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ $(LAPACK_LIBS) $(LIBS) $(LIBS) +sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp +sim_2p_incomp_reorder_LDADD = \ +$(LDADD) $(LIBS) \ +$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ +$(LAPACK_LIBS) $(LIBS) $(LIBS) + wells_example_SOURCES = wells_example.cpp if UMFPACK diff --git a/examples/sim_2p_incomp_reorder.cpp b/examples/sim_2p_incomp_reorder.cpp new file mode 100644 index 00000000..0b1ed7fc --- /dev/null +++ b/examples/sim_2p_incomp_reorder.cpp @@ -0,0 +1,220 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +{ + using namespace Opm; + + std::cout << "\n================ Test program for incompressible two-phase flow ===============\n\n"; + Opm::parameter::ParameterGroup param(argc, argv, false); + std::cout << "--------------- Reading parameters ---------------" << std::endl; + + // If we have a "deck_filename", grid and props will be read from that. + bool use_deck = param.has("deck_filename"); + boost::scoped_ptr grid; + boost::scoped_ptr props; + boost::scoped_ptr wells; + boost::scoped_ptr rock_comp; + Opm::SimulatorTimer simtimer; + Opm::TwophaseState state; + // bool check_well_controls = false; + // int max_well_control_iterations = 0; + double gravity[3] = { 0.0 }; + if (use_deck) { + std::string deck_filename = param.get("deck_filename"); + Opm::EclipseGridParser deck(deck_filename); + // Grid init + grid.reset(new Opm::GridManager(deck)); + // Rock and fluid init + const int* gc = grid->c_grid()->global_cell; + std::vector global_cell(gc, gc + grid->c_grid()->number_of_cells); + props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell)); + // Wells init. + wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + // Timer init. + if (deck.hasField("TSTEP")) { + simtimer.init(deck); + } else { + simtimer.init(param); + } + // Rock compressibility. + rock_comp.reset(new Opm::RockCompressibility(deck)); + // Gravity. + gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity; + // Init state variables (saturation and pressure). + if (param.has("init_saturation")) { + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + } else { + initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + } + } else { + // Grid init. + const int nx = param.getDefault("nx", 100); + const int ny = param.getDefault("ny", 100); + const int nz = param.getDefault("nz", 1); + const double dx = param.getDefault("dx", 1.0); + const double dy = param.getDefault("dy", 1.0); + const double dz = param.getDefault("dz", 1.0); + grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz)); + // Rock and fluid init. + props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + // Wells init. + wells.reset(new Opm::WellsManager()); + // Timer init. + simtimer.init(param); + // Rock compressibility. + rock_comp.reset(new Opm::RockCompressibility(param)); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + } + + // Warn if gravity but no density difference. + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + if (use_gravity) { + if (props->density()[0] == props->density()[1]) { + std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl; + } + } + + // Source-related variables init. + int num_cells = grid->c_grid()->number_of_cells; + std::vector totmob; + std::vector omega; // Will remain empty if no gravity. + std::vector rc; // Will remain empty if no rock compressibility. + + // Extra rock init. + std::vector porevol; + if (rock_comp->isActive()) { + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + } else { + computePorevolume(*grid->c_grid(), props->porosity(), porevol); + } + double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + + // Initialising src + std::vector src(num_cells, 0.0); + if (wells->c_wells()) { + // Do nothing, wells will be the driving force, not source terms. + // Opm::wellsToSrc(*wells->c_wells(), num_cells, src); + } else { + const double default_injection = use_gravity ? 0.0 : 0.1; + const double flow_per_sec = param.getDefault("injected_porevolumes_per_day", default_injection) + *tot_porevol_init/Opm::unit::day; + src[0] = flow_per_sec; + src[num_cells - 1] = -flow_per_sec; + } + + // Boundary conditions. + Opm::FlowBCManager bcs; + if (param.getDefault("use_pside", false)) { + int pside = param.get("pside"); + double pside_pressure = param.get("pside_pressure"); + bcs.pressureSide(*grid->c_grid(), Opm::FlowBCManager::Side(pside), pside_pressure); + } + + // Linear solver. + Opm::LinearSolverFactory linsolver(param); + + const double *grav = use_gravity ? &gravity[0] : 0; + + + Opm::SimulatorTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells->c_wells(), + src, + bcs.c_bcs(), + linsolver, + grav); + + // Warn if any parameters are unused. + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + + // Write parameters used for later reference. + // if (output) { + // param.writeParam(output_dir + "/spu_2p.param"); + // } + + WellState well_state; + well_state.init(wells->c_wells(), state); + + simulator.run(simtimer, state, well_state); +} diff --git a/opm/core/simulator/SimulatorTwophase.cpp b/opm/core/simulator/SimulatorTwophase.cpp index a83b57b7..1a617a4b 100644 --- a/opm/core/simulator/SimulatorTwophase.cpp +++ b/opm/core/simulator/SimulatorTwophase.cpp @@ -17,8 +17,40 @@ along with OPM. If not, see . */ +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + #include #include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + namespace Opm { @@ -26,87 +58,404 @@ namespace Opm class SimulatorTwophase::Impl { public: - Impl() {} - void init(const parameter::ParameterGroup& param); - void run(const SimulatorTimer& timer, const Wells& wells, TwophaseState& state, WellState& well_state); + Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity); + + void run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state); + private: -#if 0 + // Data. + + // Parameters for output. bool output_; std::string output_dir_; - int output_interval; + int output_interval_; + // Parameters for pressure solver. + int nl_pressure_maxiter_; + double nl_pressure_tolerance_; + // Parameters for transport solver. + int nl_maxiter_; + double nl_tolerance_; int num_transport_substeps_; - double gravity_[3]; - bool use_gravity_; - bool use_segregation_split = false; // - int nl_pressure_maxiter = 0; - double nl_pressure_tolerance = 0.0; - const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); - const int nl_maxiter = param.getDefault("nl_maxiter", 30); - - // boost::scoped_ptr grid_; - // boost::scoped_ptr props_; - // boost::scoped_ptr wells_; - // boost::scoped_ptr rock_comp_; - // Opm::SimulatorTimer simtimer_; - // Opm::TwophaseState state_; - // std::vector src(num_cells, 0.0); - // Opm::FlowBCManager bcs; - // Opm::LinearSolverFactory linsolver(param); + bool use_segregation_split_; + // Observed objects. + const UnstructuredGrid& grid_; + const IncompPropertiesInterface& props_; + const RockCompressibility* rock_comp_; + const Wells* wells_; + const std::vector& src_; + const FlowBoundaryConditions* bcs_; + const LinearSolverInterface& linsolver_; + const double* gravity_; + // Solvers + IncompTpfa psolver_; + TransportModelTwophase tsolver_; + // Needed by column-based gravity segregation solver. + std::vector< std::vector > columns_; + // Misc. data + std::vector allcells_; + }; + SimulatorTwophase::SimulatorTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity) + { + pimpl_.reset(new Impl(param, grid, props, rock_comp, wells, src, bcs, linsolver, gravity)); + } + + + + + + void SimulatorTwophase::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) + { + pimpl_->run(timer, state, well_state); + } + + + + static void outputState(const UnstructuredGrid& grid, + const Opm::TwophaseState& state, + const int step, + const std::string& output_dir) + { + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + if (!vtkfile) { + THROW("Failed to open " << vtkfilename.str()); + } + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + Opm::writeVtkData(grid, dm, vtkfile); + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + + + static void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir) + { + // Write water cut curve. + std::string fname = output_dir + "/watercut.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); + } + watercut.write(os); + } + + + static void outputWellReport(const Opm::WellReport& wellreport, + const std::string& output_dir) + { + // Write well report. + std::string fname = output_dir + "/wellreport.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); + } + wellreport.write(os); + } + + + + + + SimulatorTwophase::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity) + : grid_(grid), + props_(props), + rock_comp_(rock_comp), + wells_(wells), + src_(src), + bcs_(bcs), + linsolver_(linsolver), + gravity_(gravity), + psolver_(grid, props.permeability(), gravity, linsolver), + tsolver_(grid, props, 1e-9, 30) + { + // For output. + output_ = param.getDefault("output", true); + if (output_) { + output_dir_ = param.getDefault("output_dir", std::string("output")); + // Ensure that output dir exists + boost::filesystem::path fpath(output_dir_); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + output_interval_ = param.getDefault("output_interval", 1); + } + + // For pressure solver + nl_pressure_maxiter_ = param.getDefault("nl_pressure_maxiter", 10); + nl_pressure_tolerance_ = param.getDefault("nl_pressure_tolerance", 1.0); // Pascal + + // For transport solver. + nl_maxiter_ = param.getDefault("nl_maxiter", 30); + nl_tolerance_ = param.getDefault("nl_tolerance", 1e-9); + num_transport_substeps_ = param.getDefault("num_transport_substeps", 1); + use_segregation_split_ = param.getDefault("use_segregation_split", false); + if (gravity != 0 && use_segregation_split_){ + tsolver_.initGravity(gravity); + extractColumn(grid_, columns_); + } + + // Misc init. + const int num_cells = grid.number_of_cells; + allcells_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells_[cell] = cell; + } + } + + + + + void SimulatorTwophase::Impl::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) + { std::vector totmob; std::vector omega; // Will remain empty if no gravity. std::vector rc; // Will remain empty if no rock compressibility. - std::vector porevol; - double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); std::vector transport_src; - Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver); + // Initialisation. + std::vector porevol; + if (rock_comp_ && rock_comp_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); + } else { + computePorevolume(grid_, props_.porosity(), porevol); + } + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); - Opm::TransportModelTwophase tsolver(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); - typedef std::pair, std::vector > > ColMap; - ColMap columns; + // Main simulation loop. + Opm::time::StopWatch pressure_timer; + double ptime = 0.0; + Opm::time::StopWatch transport_timer; + double ttime = 0.0; + Opm::time::StopWatch total_timer; + total_timer.start(); + std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl; + double init_satvol[2] = { 0.0 }; + double satvol[2] = { 0.0 }; + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol); + std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init + << " " << init_satvol[1]/tot_porevol_init << std::endl; + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); + Opm::WellReport wellreport; + std::vector fractional_flows; + std::vector well_resflows_phase; + int num_wells = 0; + if (wells_) { + num_wells = wells_->number_of_wells; + well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); + wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates()); + } + const int num_cells = grid_.number_of_cells; + for (; !timer.done(); ++timer) { + // Report timestep and (optionally) write state to disk. + timer.report(std::cout); + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + outputState(grid_, state, timer.currentStepNum(), output_dir_); + } - std::vector allcells(num_cells); -#endif - }; + // Solve pressure. + if (gravity_) { + computeTotalMobilityOmega(props_, allcells_, state.saturation(), totmob, omega); + } else { + computeTotalMobility(props_, allcells_, state.saturation(), totmob); + } + std::vector wdp; + if (wells_) { + Opm::computeWDP(*wells_, grid_, state.saturation(), props_.density(), + gravity_ ? gravity_[2] : 0.0, true, wdp); + } + do { + pressure_timer.start(); + if (rock_comp_ && rock_comp_->isActive()) { + rc.resize(num_cells); + std::vector initial_pressure = state.pressure(); + std::vector initial_porevolume(num_cells); + computePorevolume(grid_, props_.porosity(), *rock_comp_, initial_pressure, initial_porevolume); + std::vector pressure_increment(num_cells + num_wells); + std::vector prev_pressure(num_cells + num_wells); + for (int iter = 0; iter < nl_pressure_maxiter_; ++iter) { + for (int cell = 0; cell < num_cells; ++cell) { + rc[cell] = rock_comp_->rockComp(state.pressure()[cell]); + } + computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); + std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin()); + std::copy(well_state.bhp().begin(), well_state.bhp().end(), prev_pressure.begin() + num_cells); + // prev_pressure = state.pressure(); - SimulatorTwophase::SimulatorTwophase() - { - pimpl_.reset(new Impl); - } + // compute pressure increment + psolver_.solveIncrement(totmob, omega, src_, wdp, bcs_, porevol, rc, + prev_pressure, initial_porevolume, timer.currentStepLength(), + pressure_increment); + + double max_change = 0.0; + for (int cell = 0; cell < num_cells; ++cell) { + state.pressure()[cell] += pressure_increment[cell]; + max_change = std::max(max_change, std::fabs(pressure_increment[cell])); + } + for (int well = 0; well < num_wells; ++well) { + well_state.bhp()[well] += pressure_increment[num_cells + well]; + max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well])); + } + + std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; + if (max_change < nl_pressure_tolerance_) { + break; + } + } + psolver_.computeFaceFlux(totmob, omega, src_, wdp, bcs_, state.pressure(), state.faceflux(), + well_state.bhp(), well_state.perfRates()); + } else { + psolver_.solve(totmob, omega, src_, wdp, bcs_, state.pressure(), state.faceflux(), + well_state.bhp(), well_state.perfRates()); + } + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + } while (false); + + // Process transport sources (to include bdy terms and well flows). + Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0, + wells_, well_state.perfRates(), transport_src); + + // Solve transport. + transport_timer.start(); + double stepsize = timer.currentStepLength(); + if (num_transport_substeps_ != 1) { + stepsize /= double(num_transport_substeps_); + std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; + } + for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { + tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], + stepsize, state.saturation()); + Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, injected, produced); + if (use_segregation_split_) { + tsolver_.solveGravity(columns_, &porevol[0], stepsize, state.saturation()); + } + } + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + std::cout << "Transport solver took: " << tt << " seconds." << std::endl; + ttime += tt; + + // Report volume balances. + Opm::computeSaturatedVol(porevol, state.saturation(), satvol); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + std::cout.precision(5); + const int width = 18; + std::cout << "\nVolume balance report (all numbers relative to total pore volume).\n"; + std::cout << " Saturated volumes: " + << std::setw(width) << satvol[0]/tot_porevol_init + << std::setw(width) << satvol[1]/tot_porevol_init << std::endl; + std::cout << " Injected volumes: " + << std::setw(width) << injected[0]/tot_porevol_init + << std::setw(width) << injected[1]/tot_porevol_init << std::endl; + std::cout << " Produced volumes: " + << std::setw(width) << produced[0]/tot_porevol_init + << std::setw(width) << produced[1]/tot_porevol_init << std::endl; + std::cout << " Total inj volumes: " + << std::setw(width) << tot_injected[0]/tot_porevol_init + << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl; + std::cout << " Total prod volumes: " + << std::setw(width) << tot_produced[0]/tot_porevol_init + << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl; + std::cout << " In-place + prod - inj: " + << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init + << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl; + std::cout << " Init - now - pr + inj: " + << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init + << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init + << std::endl; + std::cout.precision(8); + + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, state.saturation(), + timer.currentTime() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } + } + total_timer.stop(); + + std::cout << "\n\n================ End of simulation ===============\n" + << "Total time taken: " << total_timer.secsSinceStart() + << "\n Pressure time: " << ptime + << "\n Transport time: " << ttime << std::endl; + + if (output_) { + outputState(grid_, state, timer.currentStepNum(), output_dir_); + outputWaterCut(watercut, output_dir_); + if (wells_) { + outputWellReport(wellreport, output_dir_); + } + } - void SimulatorTwophase::init(const parameter::ParameterGroup& param) - { - pimpl_->init(param); - } - - - - void SimulatorTwophase::run(const SimulatorTimer& timer, const Wells& wells, TwophaseState& state, WellState& well_state) - { - pimpl_->run(timer, wells, state, well_state); - } - - - - void SimulatorTwophase::Impl::init(const parameter::ParameterGroup& /*param*/) - { - } - - - - void SimulatorTwophase::Impl::run(const SimulatorTimer& timer, const Wells& wells, TwophaseState& state, WellState& well_state) - { - (void) timer; - (void) wells; - (void) state; - (void) well_state; } } // namespace Opm diff --git a/opm/core/simulator/SimulatorTwophase.hpp b/opm/core/simulator/SimulatorTwophase.hpp index 636d4f94..053a3196 100644 --- a/opm/core/simulator/SimulatorTwophase.hpp +++ b/opm/core/simulator/SimulatorTwophase.hpp @@ -20,14 +20,19 @@ #ifndef OPM_SIMULATORTWOPHASE_HEADER_INCLUDED #define OPM_SIMULATORTWOPHASE_HEADER_INCLUDED -#include +#include +#include struct UnstructuredGrid; struct Wells; +struct FlowBoundaryConditions; namespace Opm { namespace parameter { class ParameterGroup; } + class IncompPropertiesInterface; + class RockCompressibility; + class LinearSolverInterface; class SimulatorTimer; class TwophaseState; class WellState; @@ -36,25 +41,53 @@ namespace Opm class SimulatorTwophase { public: - /// Default constructor. - SimulatorTwophase(); - - /// Initialise from parameters. - void init(const parameter::ParameterGroup& param); + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_pressure_tolerance (1.0) pressure solver nonlinear tolerance (in Pascal) + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] grid grid data structure + /// \param[in] props fluid and rock properties + /// \param[in] rock_comp if non-null, rock compressibility properties + /// \param[in] wells if non-null, wells data structure + /// \param[in] src source terms + /// \param[in] bcs boundary conditions, treat as all noflow if null + /// \param[in] linsolver linear solver + /// \param[in] gravity if non-null, gravity vector + SimulatorTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity); /// Run the simulation. - /// \param[in] timer governs the requested reporting timesteps - /// \param[in] wells data structure for wells - /// \param[out] state state of reservoir: pressure, fluxes - /// \param[out] well_state state of wells: bhp, perforation rates - void run(const SimulatorTimer& timer, - const Wells& wells, + /// This will run succesive timesteps until timer.done() is true. It will + /// modify the reservoir and well states. + /// \param[in,out] timer governs the requested reporting timesteps + /// \param[in,out] state state of reservoir: pressure, fluxes + /// \param[in,out] well_state state of wells: bhp, perforation rates + void run(SimulatorTimer& timer, TwophaseState& state, WellState& well_state); private: class Impl; - boost::scoped_ptr pimpl_; + // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. + boost::shared_ptr pimpl_; }; } // namespace Opm From 49641f935df0b1932ea1f1a4c20daf43f8d8afdf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 14:07:03 +0200 Subject: [PATCH 20/23] Refactor full "LDADD" statement out to single make variable. Use where appropriate. --- examples/Makefile.am | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/examples/Makefile.am b/examples/Makefile.am index 7a178da6..8c1c65c2 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -4,6 +4,13 @@ $(BOOST_CPPFLAGS) LDADD = $(top_builddir)/libopmcore.la +FULL_LDADD = \ +$(LDADD) $(LIBS) \ +$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ +$(LAPACK_LIBS) $(LIBS) $(LIBS) + + + noinst_PROGRAMS = \ refine_wells \ scaneclipsedeck \ @@ -14,16 +21,10 @@ sim_2p_incomp_reorder refine_wells_SOURCES = refine_wells.cpp sim_wateroil_SOURCES = sim_wateroil.cpp -sim_wateroil_LDADD = \ -$(LDADD) $(LIBS) \ -$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ -$(LAPACK_LIBS) $(LIBS) $(LIBS) +sim_wateroil_LDADD = $(FULL_LDADD) sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp -sim_2p_incomp_reorder_LDADD = \ -$(LDADD) $(LIBS) \ -$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ -$(LAPACK_LIBS) $(LIBS) $(LIBS) +sim_2p_incomp_reorder_LDADD = $(FULL_LDADD) wells_example_SOURCES = wells_example.cpp @@ -31,8 +32,5 @@ if UMFPACK noinst_PROGRAMS += spu_2p spu_2p_SOURCES = spu_2p.cpp -spu_2p_LDADD = \ -$(LDADD) $(LIBS) \ -$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ -$(LAPACK_LIBS) $(LIBS) $(LIBS) +spu_2p_LDADD = $(FULL_LDADD) endif From 1817d5569859b9744570aae08bc1258a5eaabe4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 14:08:44 +0200 Subject: [PATCH 21/23] Sort list of examples and PSVs. --- examples/Makefile.am | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/Makefile.am b/examples/Makefile.am index 8c1c65c2..007607c2 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -14,18 +14,18 @@ $(LAPACK_LIBS) $(LIBS) $(LIBS) noinst_PROGRAMS = \ refine_wells \ scaneclipsedeck \ +sim_2p_incomp_reorder \ sim_wateroil \ -wells_example \ -sim_2p_incomp_reorder +wells_example refine_wells_SOURCES = refine_wells.cpp -sim_wateroil_SOURCES = sim_wateroil.cpp -sim_wateroil_LDADD = $(FULL_LDADD) - sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp sim_2p_incomp_reorder_LDADD = $(FULL_LDADD) +sim_wateroil_SOURCES = sim_wateroil.cpp +sim_wateroil_LDADD = $(FULL_LDADD) + wells_example_SOURCES = wells_example.cpp if UMFPACK From a097fde7c548cbfe9043cb8553fc1e0e2f0cc422 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 14:20:12 +0200 Subject: [PATCH 22/23] Visually split description into sections. Add comments in the process. No functional changes. --- examples/Makefile.am | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/examples/Makefile.am b/examples/Makefile.am index 007607c2..7f91a6d3 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -1,15 +1,23 @@ +# Build-time flags needed to form example programs AM_CPPFLAGS = \ -I$(top_srcdir) \ $(BOOST_CPPFLAGS) +# All targets link to the library LDADD = $(top_builddir)/libopmcore.la +# Convenience make variable to avoid having to list the entire +# link-time dependency for every target that requires more than +# libopmcore.la FULL_LDADD = \ $(LDADD) $(LIBS) \ $(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ $(LAPACK_LIBS) $(LIBS) $(LIBS) - +# ---------------------------------------------------------------------- +# Declare products (i.e., the example programs). +# +# Please keep the list sorted. noinst_PROGRAMS = \ refine_wells \ @@ -18,6 +26,13 @@ sim_2p_incomp_reorder \ sim_wateroil \ wells_example +# ---------------------------------------------------------------------- +# Product constituents. Must be specified for every product that's +# built from more than a single ".c" file and/or that link to anything +# more than the OPM-Core library. +# +# Please maintain sort order from "noinst_PROGRAMS". + refine_wells_SOURCES = refine_wells.cpp sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp @@ -28,6 +43,9 @@ sim_wateroil_LDADD = $(FULL_LDADD) wells_example_SOURCES = wells_example.cpp +# ---------------------------------------------------------------------- +# Optional examples, or examples that use optional add-on components. + if UMFPACK noinst_PROGRAMS += spu_2p From 15494cd88af1766e4818ee9dc21fef1d708b2756 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 6 Jun 2012 14:42:36 +0200 Subject: [PATCH 23/23] Remove extended "LDADD" variable. Links inferred from libopmcore.la . --- examples/Makefile.am | 24 +++++------------------- 1 file changed, 5 insertions(+), 19 deletions(-) diff --git a/examples/Makefile.am b/examples/Makefile.am index 7f91a6d3..a6c29e80 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -6,14 +6,6 @@ $(BOOST_CPPFLAGS) # All targets link to the library LDADD = $(top_builddir)/libopmcore.la -# Convenience make variable to avoid having to list the entire -# link-time dependency for every target that requires more than -# libopmcore.la -FULL_LDADD = \ -$(LDADD) $(LIBS) \ -$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ -$(LAPACK_LIBS) $(LIBS) $(LIBS) - # ---------------------------------------------------------------------- # Declare products (i.e., the example programs). # @@ -33,15 +25,10 @@ wells_example # # Please maintain sort order from "noinst_PROGRAMS". -refine_wells_SOURCES = refine_wells.cpp - -sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp -sim_2p_incomp_reorder_LDADD = $(FULL_LDADD) - -sim_wateroil_SOURCES = sim_wateroil.cpp -sim_wateroil_LDADD = $(FULL_LDADD) - -wells_example_SOURCES = wells_example.cpp +refine_wells_SOURCES = refine_wells.cpp +sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp +sim_wateroil_SOURCES = sim_wateroil.cpp +wells_example_SOURCES = wells_example.cpp # ---------------------------------------------------------------------- # Optional examples, or examples that use optional add-on components. @@ -49,6 +36,5 @@ wells_example_SOURCES = wells_example.cpp if UMFPACK noinst_PROGRAMS += spu_2p -spu_2p_SOURCES = spu_2p.cpp -spu_2p_LDADD = $(FULL_LDADD) +spu_2p_SOURCES = spu_2p.cpp endif