diff --git a/examples/wells_example.cpp b/examples/wells_example.cpp index 68b54ce5..2cc9bd8e 100644 --- a/examples/wells_example.cpp +++ b/examples/wells_example.cpp @@ -3,8 +3,11 @@ #include #include #include - +#include +#include +#include #include "opm/core/newwells.h" +#include "opm/core/grid.h" int main(int argc, char** argv) { @@ -12,7 +15,7 @@ int main(int argc, char** argv) { using namespace Opm; ParameterGroup parameters( argc, argv, false ); std::string file_name = parameters.getDefault("inputdeck", "data.data"); - + // Read input file EclipseGridParser parser(file_name); std::cout << "Done!" << std::endl; @@ -22,6 +25,18 @@ int main(int argc, char** argv) { // Finally handle the wells WellsManager wells(parser, *grid.c_grid(), NULL); + std::vector global_cells(grid.c_grid()->global_cell, grid.c_grid()->global_cell + grid.c_grid()->number_of_cells); + + std::cout << "ahoi" << std::endl; + double gravity[3] = {0.0, 0.0, parameters.getDefault("gravity", 0.0)}; + IncompPropertiesFromDeck incomp_properties(parser, global_cells); + std::cout << "there" << std::endl; + + LinearSolverUmfpack umfpack_solver; + std::cout << "here" << std::endl; + IncompTpfa pressure_solver(*grid.c_grid(), incomp_properties.permeability(), + gravity, umfpack_solver, wells.c_wells()); + return 0; } diff --git a/opm/core/WellCollection.cpp b/opm/core/WellCollection.cpp index e2568f55..9d1b9f27 100644 --- a/opm/core/WellCollection.cpp +++ b/opm/core/WellCollection.cpp @@ -41,11 +41,11 @@ namespace Opm } std::tr1::shared_ptr child; - for (int i = 0; i < roots_.size(); ++i) { + for (size_t i = 0; i < roots_.size(); ++i) { if (roots_[i]->name() == child_name) { child = roots_[i]; // We've found a new parent to the previously thought root, need to remove it - for(int j = i; j < roots_.size() - 1; ++j) { + for(size_t j = i; j < roots_.size() - 1; ++j) { roots_[j] = roots_[j+1]; } @@ -63,13 +63,19 @@ namespace Opm } parent_as_group->addChild(child); - + if(child->isLeafNode()) { + leaf_nodes_.push_back(child); + } + } + + const std::vector >& WellCollection::getLeafNodes() const { + return leaf_nodes_; } WellsGroupInterface* WellCollection::findNode(std::string name) { - for (int i = 0; i < roots_.size(); i++) { + for (size_t i = 0; i < roots_.size(); i++) { WellsGroupInterface* result = roots_[i]->findGroup(name); if (result) { return result; diff --git a/opm/core/WellCollection.hpp b/opm/core/WellCollection.hpp index 53ef008e..a4cd782c 100644 --- a/opm/core/WellCollection.hpp +++ b/opm/core/WellCollection.hpp @@ -38,9 +38,14 @@ namespace Opm void addChild(std::string child, std::string parent, const EclipseGridParser& deck); + const std::vector >& getLeafNodes() const; private: // To account for the possibility of a forest std::vector > roots_; + + // This will be used to traverse the bottom nodes. + std::vector > leaf_nodes_; + WellsGroupInterface* findNode(std::string name); }; diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index 4b36acf8..17fa8a4d 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -31,13 +31,17 @@ namespace Opm : WellsGroupInterface(name, prod_spec, inj_spec) { } + + bool WellsGroupInterface::isLeafNode() const { + return false; + } WellsGroupInterface* WellsGroup::findGroup(std::string name_of_node) { if (name() == name_of_node) { return this; } else { - for (int i = 0; i < children_.size(); i++) { + for (size_t i = 0; i < children_.size(); i++) { WellsGroupInterface* result = children_[i]->findGroup(name_of_node); if (result) { return result; @@ -68,6 +72,10 @@ namespace Opm } } + + bool WellNode::isLeafNode() const { + return true; + } surface_component toSurfaceComponent(std::string type) { @@ -171,7 +179,7 @@ namespace Opm bool isWell = false; if (deck.hasField("WELSPECS")) { WELSPECS wspecs = deck.getWELSPECS(); - for (int i = 0; i < wspecs.welspecs.size(); i++) { + for (size_t i = 0; i < wspecs.welspecs.size(); i++) { if (wspecs.welspecs[i].name_ == name) { isWell = true; break; @@ -185,7 +193,7 @@ namespace Opm InjectionSpecification injection_specification; if (deck.hasField("WCONINJE")) { WCONINJE wconinje = deck.getWCONINJE(); - for (int i = 0; i < wconinje.wconinje.size(); i++) { + for (size_t i = 0; i < wconinje.wconinje.size(); i++) { if (wconinje.wconinje[i].well_ == name) { WconinjeLine line = wconinje.wconinje[i]; injection_specification.BHP_limit_ = line.BHP_limit_; @@ -199,7 +207,7 @@ namespace Opm ProductionSpecification production_specification; if (deck.hasField("WCONPROD")) { WCONPROD wconprod = deck.getWCONPROD(); - for (int i = 0; i < wconprod.wconprod.size(); i++) { + for (size_t i = 0; i < wconprod.wconprod.size(); i++) { if (wconprod.wconprod[i].well_ == name) { WconprodLine line = wconprod.wconprod[i]; production_specification.BHP_limit_ = line.BHP_limit_; @@ -216,7 +224,7 @@ namespace Opm InjectionSpecification injection_specification; if (deck.hasField("GCONINJE")) { GCONINJE gconinje = deck.getGCONINJE(); - for (int i = 0; i < gconinje.gconinje.size(); i++) { + for (size_t i = 0; i < gconinje.gconinje.size(); i++) { if (gconinje.gconinje[i].group_ == name) { GconinjeLine line = gconinje.gconinje[i]; injection_specification.injector_type_ = toSurfaceComponent(line.injector_type_); @@ -229,7 +237,7 @@ namespace Opm ProductionSpecification production_specification; if (deck.hasField("GCONPROD")) { GCONPROD gconprod = deck.getGCONPROD(); - for (int i = 0; i < gconprod.gconprod.size(); i++) { + for (size_t i = 0; i < gconprod.gconprod.size(); i++) { if (gconprod.gconprod[i].group_ == name) { GconprodLine line = gconprod.gconprod[i]; production_specification.oil_max_rate_ = line.oil_max_rate_; diff --git a/opm/core/WellsGroup.hpp b/opm/core/WellsGroup.hpp index 9d5131ce..7a9ad14d 100644 --- a/opm/core/WellsGroup.hpp +++ b/opm/core/WellsGroup.hpp @@ -17,9 +17,17 @@ namespace Opm InjectionSpecification inj_spec); virtual ~WellsGroupInterface(); + /// The unique identifier for the well or well group. const std::string& name(); + + /// Production specifications for the well or well group. const ProductionSpecification& prodSpec() const; + + /// Injection specifications for the well or well group. const InjectionSpecification& injSpec() const; + + /// \returns true if the object is a leaf node (WellNode), false otherwise. + virtual bool isLeafNode() const; /// \returns the pointer to the WellsGroupInterface with the given name. NULL if /// the name is not found.a @@ -51,6 +59,8 @@ namespace Opm InjectionSpecification inj_spec); virtual WellsGroupInterface* findGroup(std::string name_of_node); + + virtual bool isLeafNode() const; }; diff --git a/opm/core/WellsManager.cpp b/opm/core/WellsManager.cpp index 36bc12b5..2d7fb5b2 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/WellsManager.cpp @@ -505,22 +505,22 @@ namespace Opm } } - WellCollection well_collection; + WellCollection wells_; if (deck.hasField("GRUPTREE")) { std::cout << "Found gruptree" << std::endl; const GRUPTREE& gruptree = deck.getGRUPTREE(); std::map::const_iterator it = gruptree.tree.begin(); for( ; it != gruptree.tree.end(); ++it) { - well_collection.addChild(it->first, it->second, deck); + wells_.addChild(it->first, it->second, deck); } } if(deck.hasField("WELSPECS")) { WELSPECS welspecs = deck.getWELSPECS(); - for(int i = 0; i < welspecs.welspecs.size(); ++i) { + for(size_t i = 0; i < welspecs.welspecs.size(); ++i) { WelspecsLine line = welspecs.welspecs[i]; - well_collection.addChild(line.name_, line.group_, deck); + wells_.addChild(line.name_, line.group_, deck); } } diff --git a/opm/core/WellsManager.hpp b/opm/core/WellsManager.hpp index e17eaeb7..36070dd6 100644 --- a/opm/core/WellsManager.hpp +++ b/opm/core/WellsManager.hpp @@ -19,7 +19,7 @@ #ifndef OPM_WELLSMANAGER_HEADER_INCLUDED #define OPM_WELLSMANAGER_HEADER_INCLUDED - +#include struct Wells; struct UnstructuredGrid; @@ -61,6 +61,8 @@ namespace Opm WellsManager& operator=(const WellsManager& other); // The managed Wells. Wells* w_; + + WellCollection well_collection_; }; diff --git a/opm/core/fluid/BlackoilPropertiesBasic.cpp b/opm/core/fluid/BlackoilPropertiesBasic.cpp index 41c8ea70..4b8dbef2 100644 --- a/opm/core/fluid/BlackoilPropertiesBasic.cpp +++ b/opm/core/fluid/BlackoilPropertiesBasic.cpp @@ -125,7 +125,7 @@ namespace Opm double B[2]; // Must be enough since component classes do not handle more than 2. pvt_.B(1, 0, 0, B); // Compute A matrix -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { double* m = A + i*np*np; std::fill(m, m + np*np, 0.0); @@ -137,66 +137,12 @@ namespace Opm // Derivative of A matrix. if (dAdp) { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { double* m = dAdp + i*np*np; std::fill(m, m + np*np, 0.0); } } -#if 0 - // This is copied from BlackoilPropertiesFromDeck. - const int np = numPhases(); - B_.resize(n*np); - R_.resize(n*np); - if (dAdp) { - dB_.resize(n*np); - dR_.resize(n*np); - pvt_.dBdp(n, p, z, &B_[0], &dB_[0]); - pvt_.dRdp(n, p, z, &R_[0], &dR_[0]); - } else { - pvt_.B(n, p, z, &B_[0]); - pvt_.R(n, p, z, &R_[0]); - } - const int* phase_pos = pvt_.phasePosition(); - bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] && - pvt_.phaseUsed()[BlackoilPhases::Vapour]; - const int o = phase_pos[BlackoilPhases::Liquid]; - const int g = phase_pos[BlackoilPhases::Vapour]; - - // Compute A matrix -#pragma omp parallel for - for (int i = 0; i < n; ++i) { - double* m = A + i*np*np; - std::fill(m, m + np*np, 0.0); - // Diagonal entries. - for (int phase = 0; phase < np; ++phase) { - m[phase + phase*np] = 1.0/B_[i*np + phase]; - } - // Off-diagonal entries. - if (oil_and_gas) { - m[o + g*np] = R_[i*np + g]/B_[i*np + g]; - m[g + o*np] = R_[i*np + o]/B_[i*np + o]; - } - } - - // Derivative of A matrix. - if (dAdp) { -#pragma omp parallel for - for (int i = 0; i < n; ++i) { - double* m = dAdp + i*np*np; - std::fill(m, m + np*np, 0.0); - // Diagonal entries. - for (int phase = 0; phase < np; ++phase) { - m[phase + phase*np] = -dB_[i*np + phase]/B_[i*np + phase]*B_[i*np + phase]; - } - // Off-diagonal entries. - if (oil_and_gas) { - m[o + g*np] = m[g + g*np]*R_[i*np + g] + dR_[i*np + g]/B_[i*np + g]; - m[g + o*np] = m[o + o*np]*R_[i*np + o] + dR_[i*np + o]/B_[i*np + o]; - } - } - } -#endif } /// \param[in] n Number of data points. @@ -211,7 +157,7 @@ namespace Opm { const int np = numPhases(); const double* sdens = pvt_.surfaceDensities(); -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { for (int phase = 0; phase < np; ++phase) { rho[np*i + phase] = 0.0; diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 19085a41..301a3974 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -131,7 +131,7 @@ namespace Opm const int g = phase_pos[BlackoilPhases::Vapour]; // Compute A matrix -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { double* m = A + i*np*np; std::fill(m, m + np*np, 0.0); @@ -148,7 +148,7 @@ namespace Opm // Derivative of A matrix. if (dAdp) { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { double* m = dAdp + i*np*np; std::fill(m, m + np*np, 0.0); @@ -177,7 +177,7 @@ namespace Opm { const int np = numPhases(); const double* sdens = pvt_.surfaceDensities(); -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { for (int phase = 0; phase < np; ++phase) { rho[np*i + phase] = 0.0; diff --git a/opm/core/fluid/PvtPropertiesBasic.cpp b/opm/core/fluid/PvtPropertiesBasic.cpp index c2efc0d5..56e20cfe 100644 --- a/opm/core/fluid/PvtPropertiesBasic.cpp +++ b/opm/core/fluid/PvtPropertiesBasic.cpp @@ -79,7 +79,7 @@ namespace Opm { const int np = numPhases(); for (int phase = 0; phase < np; ++phase) { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_mu[np*i + phase] = viscosity_[phase]; } @@ -93,7 +93,7 @@ namespace Opm { const int np = numPhases(); for (int phase = 0; phase < np; ++phase) { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[np*i + phase] = formation_volume_factor_[phase]; } @@ -108,7 +108,7 @@ namespace Opm { const int np = numPhases(); for (int phase = 0; phase < np; ++phase) { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[np*i + phase] = formation_volume_factor_[phase]; output_dBdp[np*i + phase] = 0.0; diff --git a/opm/core/fluid/RockBasic.cpp b/opm/core/fluid/RockBasic.cpp index 4bd2942d..a1423b76 100644 --- a/opm/core/fluid/RockBasic.cpp +++ b/opm/core/fluid/RockBasic.cpp @@ -42,7 +42,7 @@ namespace Opm permeability_.clear(); const int dsq = dimensions*dimensions; permeability_.resize(num_cells*dsq, 0.0); -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < num_cells; ++i) { for (int d = 0; d < dimensions; ++d) { permeability_[dsq*i + dimensions*d + d] = perm; diff --git a/opm/core/fluid/SaturationPropsBasic.cpp b/opm/core/fluid/SaturationPropsBasic.cpp index 5bbdd939..0f97e7c0 100644 --- a/opm/core/fluid/SaturationPropsBasic.cpp +++ b/opm/core/fluid/SaturationPropsBasic.cpp @@ -71,13 +71,13 @@ namespace Opm const double* s, double* kr, double* dkrds, Fun fun) { if (dkrds == 0) { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n*np; ++i) { kr[i] = fun.kr(s[i]); } return; } -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { std::fill(dkrds + i*np*np, dkrds + (i+1)*np*np, 0.0); for (int phase = 0; phase < np; ++phase) { diff --git a/opm/core/fluid/SaturationPropsFromDeck.cpp b/opm/core/fluid/SaturationPropsFromDeck.cpp index e70cf937..2075b999 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.cpp +++ b/opm/core/fluid/SaturationPropsFromDeck.cpp @@ -111,12 +111,12 @@ namespace Opm { const int np = phase_usage_.num_phases; if (dkrds) { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i); } } else { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { evalKr(s + np*i, kr + np*i); } @@ -142,12 +142,12 @@ namespace Opm { const int np = phase_usage_.num_phases; if (dpcds) { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i); } } else { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { evalPc(s + np*i, pc + np*i); } diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp index 9b057040..9fdaf445 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp +++ b/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp @@ -134,7 +134,7 @@ namespace Opm data1_.resize(n); for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { props_[phase]->mu(n, p, z, &data1_[0]); -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_mu[phase_usage_.num_phases*i + phase] = data1_[i]; } @@ -149,7 +149,7 @@ namespace Opm data1_.resize(n); for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { props_[phase]->B(n, p, z, &data1_[0]); -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[phase_usage_.num_phases*i + phase] = data1_[i]; } @@ -166,7 +166,7 @@ namespace Opm data2_.resize(n); for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { props_[phase]->dBdp(n, p, z, &data1_[0], &data2_[0]); -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[phase_usage_.num_phases*i + phase] = data1_[i]; output_dBdp[phase_usage_.num_phases*i + phase] = data2_[i]; @@ -183,7 +183,7 @@ namespace Opm data1_.resize(n); for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { props_[phase]->R(n, p, z, &data1_[0]); -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_R[phase_usage_.num_phases*i + phase] = data1_[i]; } @@ -200,7 +200,7 @@ namespace Opm data2_.resize(n); for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { props_[phase]->dRdp(n, p, z, &data1_[0], &data2_[0]); -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_R[phase_usage_.num_phases*i + phase] = data1_[i]; output_dRdp[phase_usage_.num_phases*i + phase] = data2_[i]; diff --git a/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp b/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp index c1ed11b8..4bf20897 100644 --- a/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp +++ b/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp @@ -81,7 +81,7 @@ namespace Opm double* output_B) const { if (comp_) { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { // Computing a polynomial approximation to the exponential. double x = comp_*(p[i] - ref_press_); @@ -100,7 +100,7 @@ namespace Opm { B(n, p, 0, output_B); if (comp_) { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_dBdp[i] = -comp_*output_B[i]; } diff --git a/opm/core/fluid/blackoil/SinglePvtDead.cpp b/opm/core/fluid/blackoil/SinglePvtDead.cpp index 431b55b3..3d1a2740 100644 --- a/opm/core/fluid/blackoil/SinglePvtDead.cpp +++ b/opm/core/fluid/blackoil/SinglePvtDead.cpp @@ -75,7 +75,7 @@ namespace Opm const double* /*z*/, double* output_mu) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_mu[i] = viscosity_(p[i]); } @@ -86,7 +86,7 @@ namespace Opm const double* /*z*/, double* output_B) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[i] = 1.0/one_over_B_(p[i]); } @@ -99,7 +99,7 @@ namespace Opm double* output_dBdp) const { B(n, p, 0, output_B); -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { double Bg = output_B[i]; output_dBdp[i] = -Bg*Bg*one_over_B_.derivative(p[i]); diff --git a/opm/core/fluid/blackoil/SinglePvtLiveGas.cpp b/opm/core/fluid/blackoil/SinglePvtLiveGas.cpp index 048bc834..2526d34c 100644 --- a/opm/core/fluid/blackoil/SinglePvtLiveGas.cpp +++ b/opm/core/fluid/blackoil/SinglePvtLiveGas.cpp @@ -91,7 +91,7 @@ namespace Opm const double* z, double* output_mu) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_mu[i] = miscible_gas(p[i], z + num_phases_*i, 2, false); } @@ -104,7 +104,7 @@ namespace Opm const double* z, double* output_B) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[i] = evalB(p[i], z + num_phases_*i); } @@ -119,7 +119,7 @@ namespace Opm double* output_B, double* output_dBdp) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { evalBDeriv(p[i], z + num_phases_*i, output_B[i], output_dBdp[i]); } @@ -132,7 +132,7 @@ namespace Opm const double* z, double* output_R) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_R[i] = evalR(p[i], z + num_phases_*i); } @@ -147,7 +147,7 @@ namespace Opm double* output_R, double* output_dRdp) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { evalRDeriv(p[i], z + num_phases_*i, output_R[i], output_dRdp[i]); } diff --git a/opm/core/fluid/blackoil/SinglePvtLiveOil.cpp b/opm/core/fluid/blackoil/SinglePvtLiveOil.cpp index 2d42997d..2a393839 100644 --- a/opm/core/fluid/blackoil/SinglePvtLiveOil.cpp +++ b/opm/core/fluid/blackoil/SinglePvtLiveOil.cpp @@ -168,7 +168,7 @@ namespace Opm const double* z, double* output_mu) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_mu[i] = miscible_oil(p[i], z + num_phases_*i, 2, false); } @@ -181,7 +181,7 @@ namespace Opm const double* z, double* output_B) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_B[i] = evalB(p[i], z + num_phases_*i); } @@ -196,7 +196,7 @@ namespace Opm double* output_B, double* output_dBdp) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { evalBDeriv(p[i], z + num_phases_*i, output_B[i], output_dBdp[i]); } @@ -209,7 +209,7 @@ namespace Opm const double* z, double* output_R) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { output_R[i] = evalR(p[i], z + num_phases_*i); } @@ -224,7 +224,7 @@ namespace Opm double* output_R, double* output_dRdp) const { -#pragma omp parallel for +// #pragma omp parallel for for (int i = 0; i < n; ++i) { evalRDeriv(p[i], z + num_phases_*i, output_R[i], output_dRdp[i]); } diff --git a/opm/core/pressure/IncompTpfa.cpp b/opm/core/pressure/IncompTpfa.cpp index 9417b7eb..8ec414d2 100644 --- a/opm/core/pressure/IncompTpfa.cpp +++ b/opm/core/pressure/IncompTpfa.cpp @@ -38,14 +38,18 @@ namespace Opm /// and N == g.number_of_cells. /// \param[in] gravity Gravity vector. If nonzero, the array should /// have D elements. + /// \param[in] wells The wells argument. Will be used in solution, + /// is ignored if NULL IncompTpfa::IncompTpfa(const UnstructuredGrid& g, const double* permeability, const double* gravity, - const LinearSolverInterface& linsolver) + const LinearSolverInterface& linsolver, + const struct Wells* wells) : grid_(g), linsolver_(linsolver), htrans_(g.cell_facepos[ g.number_of_cells ]), - trans_ (g.number_of_faces) + trans_ (g.number_of_faces), + wells_(wells) { UnstructuredGrid* gg = const_cast(&grid_); tpfa_htrans_compute(gg, permeability, &htrans_[0]); @@ -110,7 +114,7 @@ namespace Opm } } - ifs_tpfa_forces F = { NULL, NULL, NULL, NULL, NULL }; + ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL }; if (! src.empty()) { F.src = &src[0]; } F.bc = bcs; @@ -172,7 +176,7 @@ namespace Opm } } - ifs_tpfa_forces F = { NULL, NULL, NULL, NULL, NULL }; + ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL }; if (! src.empty()) { F.src = &src[0]; } F.bc = bcs; diff --git a/opm/core/pressure/IncompTpfa.hpp b/opm/core/pressure/IncompTpfa.hpp index 58b53383..01eb6627 100644 --- a/opm/core/pressure/IncompTpfa.hpp +++ b/opm/core/pressure/IncompTpfa.hpp @@ -25,6 +25,7 @@ struct UnstructuredGrid; struct ifs_tpfa_data; +struct Wells; struct FlowBoundaryConditions; namespace Opm @@ -47,10 +48,13 @@ namespace Opm /// and N == g.number_of_cells. /// \param[in] gravity Gravity vector. If nonzero, the array should /// have D elements. + /// \param[in] wells The wells argument. Will be used in solution, + /// is ignored if NULL IncompTpfa(const UnstructuredGrid& g, const double* permeability, const double* gravity, - const LinearSolverInterface& linsolver); + const LinearSolverInterface& linsolver, + const struct Wells* wells = 0); /// Destructor. ~IncompTpfa(); @@ -113,7 +117,8 @@ namespace Opm ::std::vector trans_ ; ::std::vector gpress_; ::std::vector gpress_omegaweighted_; - + + const struct Wells* wells_; struct ifs_tpfa_data* h_; }; diff --git a/opm/core/utility/writeVtkData.cpp b/opm/core/utility/writeVtkData.cpp index 3ec8fbae..59b4d8dc 100644 --- a/opm/core/utility/writeVtkData.cpp +++ b/opm/core/utility/writeVtkData.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include