From 3e449a31d6f3253dda30e970a255fc1bc44fb180 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 9 Jun 2015 11:25:36 +0200 Subject: [PATCH 01/15] adjusting some spaces related to brackets. No functions change. --- opm/polymer/PolymerProperties.hpp | 6 ++-- .../fullyimplicit/BlackoilPolymerModel.hpp | 4 +-- .../BlackoilPolymerModel_impl.hpp | 28 ++++++++++--------- 3 files changed, 20 insertions(+), 18 deletions(-) diff --git a/opm/polymer/PolymerProperties.hpp b/opm/polymer/PolymerProperties.hpp index 70e2867f6..723bb5fde 100644 --- a/opm/polymer/PolymerProperties.hpp +++ b/opm/polymer/PolymerProperties.hpp @@ -156,7 +156,7 @@ namespace Opm plyshlog_ = deck->hasKeyword("PLYSHLOG"); - if(plyshlog_){ + if (plyshlog_) { // Assuming NTPVT == 1 always due to the limitation of the parser const auto& plyshlogTable = eclipseState->getPlyshlogTables()[0]; @@ -174,14 +174,14 @@ namespace Opm plyshlog_ref_conc_ = plyshlogTable.getRefPolymerConcentration(); - if(plyshlogTable.hasRefSalinity()) { + if (plyshlogTable.hasRefSalinity()) { has_plyshlog_ref_salinity_ = true; plyshlog_ref_salinity_ = plyshlogTable.getRefSalinity(); } else { has_plyshlog_ref_salinity_ = false; } - if(plyshlogTable.hasRefTemperature()) { + if (plyshlogTable.hasRefTemperature()) { has_plyshlog_ref_temp_ = true; plyshlog_ref_temp_ = plyshlogTable.getRefTemperature(); } else { diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index cc8be471b..aa6784554 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -293,10 +293,10 @@ namespace Opm { /// Finding the intersection point of a line segment and a line. /// return true, if found. - bool findIntersection (Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point); + bool findIntersection(Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point); /// Computing the shear multiplier based on the water velocity/shear rate with PLYSHLOG keyword - bool computeShearMultLog( std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult); + bool computeShearMultLog(std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult); /// Computing the water velocity without shear-thinning for the cell faces. /// The water velocity will be used for shear-thinning calculation. diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index a7d16a532..c028be244 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -895,7 +895,7 @@ namespace Opm { template bool - BlackoilPolymerModel::findIntersection (Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point) + BlackoilPolymerModel::findIntersection(Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point) { const double x1 = line_segment1[0].x; @@ -910,12 +910,14 @@ namespace Opm { const double d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4); - if (d == 0.) { return false; } + if (d == 0.) { + return false; + } const double x = ((x3 - x4) * (x1 * y2 - y1 * x2) - (x1 - x2) * (x3 * y4 - y3 * x4)) / d; const double y = ((y3 - y4) * (x1 * y2 - y1 * x2) - (y1 - y2) * (x3 * y4 - y3 * x4)) / d; - if( x >= std::min(x1,x2) && x <= std::max(x1,x2) ){ + if (x >= std::min(x1,x2) && x <= std::max(x1,x2)) { intersection_point.x = x; intersection_point.y = y; return true; @@ -926,7 +928,7 @@ namespace Opm { template bool - BlackoilPolymerModel::computeShearMultLog( std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult) + BlackoilPolymerModel::computeShearMultLog(std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult) { double refConcentration = polymer_props_ad_.plyshlogRefConc(); @@ -942,7 +944,7 @@ namespace Opm { logShearVRF.resize(shear_water_vel.size()); // converting the table using the reference condition - for(int i = 0; i < shear_vrf.size(); ++i){ + for (int i = 0; i < shear_vrf.size(); ++i) { shear_vrf[i] = (refViscMult * shear_vrf[i] - 1.) / (refViscMult - 1); logShearWaterVel[i] = std::log(shear_water_vel[i]); } @@ -954,14 +956,14 @@ namespace Opm { const double maxShearVel = shear_water_vel.back(); const double epsilon = std::sqrt(std::numeric_limits::epsilon()); - for(int i = 0; i < water_vel.size(); ++i){ + for (int i = 0; i < water_vel.size(); ++i) { - if( visc_mult[i] - 1. < epsilon || std::abs(water_vel[i]) < minShearVel ) { + if (visc_mult[i] - 1. < epsilon || std::abs(water_vel[i]) < minShearVel) { shear_mult[i] = 1.0; continue; } - for(int j = 0; j < shear_vrf.size(); ++j){ + for (int j = 0; j < shear_vrf.size(); ++j) { logShearVRF[j] = (1 + (visc_mult[i] - 1.0) * shear_vrf[j]) / visc_mult[i]; logShearVRF[j] = std::log(logShearVRF[j]); } @@ -972,7 +974,7 @@ namespace Opm { int iIntersection; // finding the intersection on the iIntersectionth table segment bool foundSegment = false; - for(iIntersection = 0; iIntersection < shear_vrf.size() - 1; ++iIntersection){ + for (iIntersection = 0; iIntersection < shear_vrf.size() - 1; ++iIntersection) { double temp1 = logShearVRF[iIntersection] + logShearWaterVel[iIntersection] - logWaterVelO; double temp2 = logShearVRF[iIntersection + 1] + logShearWaterVel[iIntersection + 1] - logWaterVelO; @@ -985,7 +987,7 @@ namespace Opm { } } - if(foundSegment == true){ + if (foundSegment == true) { Point2D lineSegment[2]; lineSegment[0] = Point2D{logShearWaterVel[iIntersection], logShearVRF[iIntersection]}; lineSegment[1] = Point2D{logShearWaterVel[iIntersection + 1], logShearVRF[iIntersection + 1]}; @@ -998,13 +1000,13 @@ namespace Opm { bool foundIntersection = findIntersection(lineSegment, line, intersectionPoint); - if(foundIntersection){ + if (foundIntersection) { shear_mult[i] = std::exp(intersectionPoint.y); - }else{ + } else { std::cerr << " failed in finding the solution for shear-thinning multiplier " << std::endl; return false; // failed in finding the solution. } - }else{ + } else { if (std::abs(water_vel[i]) < maxShearVel) { std::cout << " the veclocity is " << water_vel[i] << std::endl; std::cout << " max shear velocity is " << maxShearVel << std::endl; From e44ef196acba13bf076f7c297b0b1e3d666316b4 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 10 Jun 2015 12:51:58 +0200 Subject: [PATCH 02/15] moving the intersection calculation to seperate point2D class. under the namespae Opm::detail --- opm/polymer/Point2D.hpp | 102 ++++++++++++++++++ .../fullyimplicit/BlackoilPolymerModel.hpp | 10 -- .../BlackoilPolymerModel_impl.hpp | 53 ++------- 3 files changed, 112 insertions(+), 53 deletions(-) create mode 100644 opm/polymer/Point2D.hpp diff --git a/opm/polymer/Point2D.hpp b/opm/polymer/Point2D.hpp new file mode 100644 index 000000000..49549bef5 --- /dev/null +++ b/opm/polymer/Point2D.hpp @@ -0,0 +1,102 @@ +/* + Copyright 2015 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_POINT2DINTERSECTION_INCLUDED +#define OPM_POINT2DINTERSECTION_INCLUDED + +namespace Opm { + + namespace detail { + class Point2D + { + public: + Point2D(const double xi, const double yi) + : x_(xi), + y_(yi) { + + } + + Point2D() + { + } + + const double getX() const + { + return x_; + } + + const double getY() const + { + return y_; + } + + void setX(const double x) + { + x_ = x; + } + + void setY(const double y){ + y_ = y; + } + + /// Finding the intersection point of a line segment and a line. + /// return true, if found. + static bool findIntersection(Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point) + { + + const double x1 = line_segment1[0].getX(); + const double y1 = line_segment1[0].getY(); + const double x2 = line_segment1[1].getX(); + const double y2 = line_segment1[1].getY(); + + const double x3 = line2[0].getX(); + const double y3 = line2[0].getY(); + const double x4 = line2[1].getX(); + const double y4 = line2[1].getY(); + + const double d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4); + + if (d == 0.) { + return false; + } + + const double x = ((x3 - x4) * (x1 * y2 - y1 * x2) - (x1 - x2) * (x3 * y4 - y3 * x4)) / d; + const double y = ((y3 - y4) * (x1 * y2 - y1 * x2) - (y1 - y2) * (x3 * y4 - y3 * x4)) / d; + + if (x >= std::min(x1,x2) && x <= std::max(x1,x2)) { + intersection_point.setX(x); + intersection_point.setY(y); + return true; + } else { + return false; + } + } + + private: + double x_; + double y_; + + }; // class Point2D + + } // namespace detail + +} // namespace Opm + +#endif // OPM_POINT2DINTERSECTION_INCLUDED + diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index aa6784554..1b8af9773 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -285,16 +285,6 @@ namespace Opm { int nc, int nw) const; - - struct Point2D { - double x; - double y; - }; - - /// Finding the intersection point of a line segment and a line. - /// return true, if found. - bool findIntersection(Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point); - /// Computing the shear multiplier based on the water velocity/shear rate with PLYSHLOG keyword bool computeShearMultLog(std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult); diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index c028be244..05df70c93 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -25,6 +25,7 @@ #define OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED #include +#include #include #include @@ -892,40 +893,6 @@ namespace Opm { extraAddWellEq(state, xw, cq_ps, cmix_s, cqt_is, well_cells); } - - template - bool - BlackoilPolymerModel::findIntersection(Point2D line_segment1[2], Point2D line2[2], Point2D& intersection_point) - { - - const double x1 = line_segment1[0].x; - const double y1 = line_segment1[0].y; - const double x2 = line_segment1[1].x; - const double y2 = line_segment1[1].y; - - const double x3 = line2[0].x; - const double y3 = line2[0].y; - const double x4 = line2[1].x; - const double y4 = line2[1].y; - - const double d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4); - - if (d == 0.) { - return false; - } - - const double x = ((x3 - x4) * (x1 * y2 - y1 * x2) - (x1 - x2) * (x3 * y4 - y3 * x4)) / d; - const double y = ((y3 - y4) * (x1 * y2 - y1 * x2) - (y1 - y2) * (x3 * y4 - y3 * x4)) / d; - - if (x >= std::min(x1,x2) && x <= std::max(x1,x2)) { - intersection_point.x = x; - intersection_point.y = y; - return true; - } else { - return false; - } - } - template bool BlackoilPolymerModel::computeShearMultLog(std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult) @@ -988,20 +955,20 @@ namespace Opm { } if (foundSegment == true) { - Point2D lineSegment[2]; - lineSegment[0] = Point2D{logShearWaterVel[iIntersection], logShearVRF[iIntersection]}; - lineSegment[1] = Point2D{logShearWaterVel[iIntersection + 1], logShearVRF[iIntersection + 1]}; + detail::Point2D lineSegment[2]; + lineSegment[0] = detail::Point2D{logShearWaterVel[iIntersection], logShearVRF[iIntersection]}; + lineSegment[1] = detail::Point2D{logShearWaterVel[iIntersection + 1], logShearVRF[iIntersection + 1]}; - Point2D line[2]; - line[0] = Point2D{0, logWaterVelO}; - line[1] = Point2D{logWaterVelO, 0}; + detail::Point2D line[2]; + line[0] = detail::Point2D{0, logWaterVelO}; + line[1] = detail::Point2D{logWaterVelO, 0}; - Point2D intersectionPoint; + detail::Point2D intersectionPoint; - bool foundIntersection = findIntersection(lineSegment, line, intersectionPoint); + bool foundIntersection = detail::Point2D::findIntersection(lineSegment, line, intersectionPoint); if (foundIntersection) { - shear_mult[i] = std::exp(intersectionPoint.y); + shear_mult[i] = std::exp(intersectionPoint.getY()); } else { std::cerr << " failed in finding the solution for shear-thinning multiplier " << std::endl; return false; // failed in finding the solution. From 8d47acd6bfbb5b3beb4bb6e23edd81064800dd39 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 10 Jun 2015 13:22:43 +0200 Subject: [PATCH 03/15] Satisfy interface requirements of SimulatorBase. This is done with dummy types since the class in question (SimulatorFullyImplicitCompressiblePolymer) does not use a Model/Solver split. --- .../fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp | 3 +++ .../SimulatorFullyImplicitCompressiblePolymer.hpp | 5 +++++ 2 files changed, 8 insertions(+) diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp index e86b34769..293ed06f1 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp @@ -88,6 +88,9 @@ namespace Opm { int newtonIterations() const; int linearIterations() const; + /// Not used by this class except to satisfy interface requirements. + typedef parameter::ParameterGroup SolverParameters; + private: typedef AutoDiffBlock ADB; typedef ADB::V V; diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp index e3abe9654..c6a9a9110 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp @@ -79,6 +79,11 @@ namespace Opm typedef BlackoilOutputWriter OutputWriter; typedef GridT Grid; typedef FullyImplicitCompressiblePolymerSolver Solver; + /// Dummy class, this Solver does not use a Model. + struct Model + { + typedef parameter::ParameterGroup ModelParameters; + }; }; /// Class collecting all necessary components for a two-phase simulation. From 81d9fe7a554d3b0f1289a4b0b5d94c96d58595ff Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 10 Jun 2015 13:36:37 +0200 Subject: [PATCH 04/15] moving function computeShearMultLog to class PolymerProperties --- opm/polymer/PolymerProperties.cpp | 95 ++++++++++++++++ opm/polymer/PolymerProperties.hpp | 4 + .../fullyimplicit/BlackoilPolymerModel.hpp | 3 - .../BlackoilPolymerModel_impl.hpp | 101 +----------------- opm/polymer/fullyimplicit/PolymerPropsAd.cpp | 8 ++ opm/polymer/fullyimplicit/PolymerPropsAd.hpp | 8 ++ 6 files changed, 117 insertions(+), 102 deletions(-) diff --git a/opm/polymer/PolymerProperties.cpp b/opm/polymer/PolymerProperties.cpp index bae878c32..87b8e8b11 100644 --- a/opm/polymer/PolymerProperties.cpp +++ b/opm/polymer/PolymerProperties.cpp @@ -20,6 +20,7 @@ #include #include +#include #include #include #include @@ -408,4 +409,98 @@ namespace Opm dmc_dc = 0.; } } + + bool PolymerProperties::computeShearMultLog(std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult) const + { + + double refConcentration = plyshlogRefConc(); + double refViscMult = viscMult(refConcentration); + + std::vector shear_water_vel = shearWaterVelocity(); + std::vector shear_vrf = shearViscosityReductionFactor(); + + std::vector logShearWaterVel; + std::vector logShearVRF; + + logShearWaterVel.resize(shear_water_vel.size()); + logShearVRF.resize(shear_water_vel.size()); + + // converting the table using the reference condition + for (size_t i = 0; i < shear_vrf.size(); ++i) { + shear_vrf[i] = (refViscMult * shear_vrf[i] - 1.) / (refViscMult - 1); + logShearWaterVel[i] = std::log(shear_water_vel[i]); + } + + shear_mult.resize(water_vel.size()); + + // the mimum velocity to apply the shear-thinning + const double minShearVel = shear_water_vel[0]; + const double maxShearVel = shear_water_vel.back(); + const double epsilon = std::sqrt(std::numeric_limits::epsilon()); + + for (size_t i = 0; i < water_vel.size(); ++i) { + + if (visc_mult[i] - 1. < epsilon || std::abs(water_vel[i]) < minShearVel) { + shear_mult[i] = 1.0; + continue; + } + + for (size_t j = 0; j < shear_vrf.size(); ++j) { + logShearVRF[j] = (1 + (visc_mult[i] - 1.0) * shear_vrf[j]) / visc_mult[i]; + logShearVRF[j] = std::log(logShearVRF[j]); + } + + // const double logWaterVelO = std::log(water_vel[i]); + const double logWaterVelO = std::log(std::abs(water_vel[i])); + + size_t iIntersection; // finding the intersection on the iIntersectionth table segment + bool foundSegment = false; + + for (iIntersection = 0; iIntersection < shear_vrf.size() - 1; ++iIntersection) { + + double temp1 = logShearVRF[iIntersection] + logShearWaterVel[iIntersection] - logWaterVelO; + double temp2 = logShearVRF[iIntersection + 1] + logShearWaterVel[iIntersection + 1] - logWaterVelO; + + // ignore the cases the temp1 or temp2 is zero first for simplicity. + // several more complicated cases remain to be implemented. + if( temp1 * temp2 < 0.){ + foundSegment = true; + break; + } + } + + if (foundSegment == true) { + detail::Point2D lineSegment[2]; + lineSegment[0] = detail::Point2D{logShearWaterVel[iIntersection], logShearVRF[iIntersection]}; + lineSegment[1] = detail::Point2D{logShearWaterVel[iIntersection + 1], logShearVRF[iIntersection + 1]}; + + detail::Point2D line[2]; + line[0] = detail::Point2D{0, logWaterVelO}; + line[1] = detail::Point2D{logWaterVelO, 0}; + + detail::Point2D intersectionPoint; + + bool foundIntersection = detail::Point2D::findIntersection(lineSegment, line, intersectionPoint); + + if (foundIntersection) { + shear_mult[i] = std::exp(intersectionPoint.getY()); + } else { + std::cerr << " failed in finding the solution for shear-thinning multiplier " << std::endl; + return false; // failed in finding the solution. + } + } else { + if (std::abs(water_vel[i]) < maxShearVel) { + std::cout << " the veclocity is " << water_vel[i] << std::endl; + std::cout << " max shear velocity is " << maxShearVel << std::endl; + std::cerr << " something wrong happend in finding segment" << std::endl; + return false; + } else { + shear_mult[i] = std::exp(logShearVRF.back()); + } + } + + } + + return true; + } } diff --git a/opm/polymer/PolymerProperties.hpp b/opm/polymer/PolymerProperties.hpp index 723bb5fde..7feff246c 100644 --- a/opm/polymer/PolymerProperties.hpp +++ b/opm/polymer/PolymerProperties.hpp @@ -326,6 +326,9 @@ namespace Opm void computeMcBoth(const double& c, double& mc, double& dmc_dc, bool if_with_der) const; + /// Computing the shear multiplier based on the water velocity/shear rate with PLYSHLOG keyword + bool computeShearMultLog(std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult) const; + private: double c_max_; double mix_param_; @@ -366,6 +369,7 @@ namespace Opm double& deff_relperm_wat_ds, double& deff_relperm_wat_dc, bool if_with_der) const; + }; } // namespace Opm diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 1b8af9773..c855707f5 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -285,9 +285,6 @@ namespace Opm { int nc, int nw) const; - /// Computing the shear multiplier based on the water velocity/shear rate with PLYSHLOG keyword - bool computeShearMultLog(std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult); - /// Computing the water velocity without shear-thinning for the cell faces. /// The water velocity will be used for shear-thinning calculation. void computeWaterShearVelocityFaces(const V& transi, const std::vector& kr, diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 05df70c93..80ee16d67 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -25,7 +25,6 @@ #define OPM_BLACKOILPOLYMERMODEL_IMPL_HEADER_INCLUDED #include -#include #include #include @@ -288,7 +287,7 @@ namespace Opm { std::vector visc_mult; computeWaterShearVelocityFaces(transi, kr, state.canonical_phase_pressures, state, water_vel, visc_mult); - if(!computeShearMultLog(water_vel, visc_mult, shear_mult_faces_)) { + if(!polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_)) { // std::cerr << " failed in calculating the shear-multiplier " << std::endl; OPM_THROW(std::runtime_error, " failed in calculating the shear-multiplier. "); } @@ -578,7 +577,7 @@ namespace Opm { computeWaterShearVelocityWells(state, well_state, aliveWells, water_vel_wells, visc_mult_wells); - if (!computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_)) { + if (!polymer_props_ad_.computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_)) { // std::cout << " failed in calculating the shear factors for wells " << std::endl; OPM_THROW(std::runtime_error, " failed in calculating the shear factors for wells "); } @@ -893,102 +892,6 @@ namespace Opm { extraAddWellEq(state, xw, cq_ps, cmix_s, cqt_is, well_cells); } - template - bool - BlackoilPolymerModel::computeShearMultLog(std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult) - { - - double refConcentration = polymer_props_ad_.plyshlogRefConc(); - double refViscMult = polymer_props_ad_.viscMult(refConcentration); - - std::vector shear_water_vel = polymer_props_ad_.shearWaterVelocity(); - std::vector shear_vrf = polymer_props_ad_.shearViscosityReductionFactor(); - - std::vector logShearWaterVel; - std::vector logShearVRF; - - logShearWaterVel.resize(shear_water_vel.size()); - logShearVRF.resize(shear_water_vel.size()); - - // converting the table using the reference condition - for (int i = 0; i < shear_vrf.size(); ++i) { - shear_vrf[i] = (refViscMult * shear_vrf[i] - 1.) / (refViscMult - 1); - logShearWaterVel[i] = std::log(shear_water_vel[i]); - } - - shear_mult.resize(water_vel.size()); - - // the mimum velocity to apply the shear-thinning - const double minShearVel = shear_water_vel[0]; - const double maxShearVel = shear_water_vel.back(); - const double epsilon = std::sqrt(std::numeric_limits::epsilon()); - - for (int i = 0; i < water_vel.size(); ++i) { - - if (visc_mult[i] - 1. < epsilon || std::abs(water_vel[i]) < minShearVel) { - shear_mult[i] = 1.0; - continue; - } - - for (int j = 0; j < shear_vrf.size(); ++j) { - logShearVRF[j] = (1 + (visc_mult[i] - 1.0) * shear_vrf[j]) / visc_mult[i]; - logShearVRF[j] = std::log(logShearVRF[j]); - } - - // const double logWaterVelO = std::log(water_vel[i]); - const double logWaterVelO = std::log(std::abs(water_vel[i])); - - int iIntersection; // finding the intersection on the iIntersectionth table segment - bool foundSegment = false; - - for (iIntersection = 0; iIntersection < shear_vrf.size() - 1; ++iIntersection) { - - double temp1 = logShearVRF[iIntersection] + logShearWaterVel[iIntersection] - logWaterVelO; - double temp2 = logShearVRF[iIntersection + 1] + logShearWaterVel[iIntersection + 1] - logWaterVelO; - - // ignore the cases the temp1 or temp2 is zero first for simplicity. - // several more complicated cases remain to be implemented. - if( temp1 * temp2 < 0.){ - foundSegment = true; - break; - } - } - - if (foundSegment == true) { - detail::Point2D lineSegment[2]; - lineSegment[0] = detail::Point2D{logShearWaterVel[iIntersection], logShearVRF[iIntersection]}; - lineSegment[1] = detail::Point2D{logShearWaterVel[iIntersection + 1], logShearVRF[iIntersection + 1]}; - - detail::Point2D line[2]; - line[0] = detail::Point2D{0, logWaterVelO}; - line[1] = detail::Point2D{logWaterVelO, 0}; - - detail::Point2D intersectionPoint; - - bool foundIntersection = detail::Point2D::findIntersection(lineSegment, line, intersectionPoint); - - if (foundIntersection) { - shear_mult[i] = std::exp(intersectionPoint.getY()); - } else { - std::cerr << " failed in finding the solution for shear-thinning multiplier " << std::endl; - return false; // failed in finding the solution. - } - } else { - if (std::abs(water_vel[i]) < maxShearVel) { - std::cout << " the veclocity is " << water_vel[i] << std::endl; - std::cout << " max shear velocity is " << maxShearVel << std::endl; - std::cerr << " something wrong happend in finding segment" << std::endl; - return false; - } else { - shear_mult[i] = std::exp(logShearVRF.back()); - } - } - - } - - return true; - } - template void BlackoilPolymerModel::computeWaterShearVelocityFaces(const V& transi, const std::vector& kr, diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp index 0fb23f6e7..37dc56ba4 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.cpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.cpp @@ -317,4 +317,12 @@ namespace Opm { return krw / rk; } + + bool + PolymerPropsAd::computeShearMultLog(std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult) const + { + return polymer_props_.computeShearMultLog(water_vel, visc_mult, shear_mult); + } + + }// namespace Opm diff --git a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp index 07ddc7a11..7bf90ea3c 100644 --- a/opm/polymer/fullyimplicit/PolymerPropsAd.hpp +++ b/opm/polymer/fullyimplicit/PolymerPropsAd.hpp @@ -130,6 +130,14 @@ namespace Opm { ADB effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw) const; + /// \param[in] water_vel Array of the n values of water velocity or shear rate. + /// \param[in] visc_mult Array of the n values of the viscosity multiplier from PLYVISC table. + /// \parma[out] shear_mult Array of the n values of calculated shear multiplier with PLYSHLOG keyword. + /// \return TRUE if the calculation of shear multiplier is sucessful, + /// FALSE if the calculation of shear multplier is failed. + bool computeShearMultLog(std::vector& water_vel, std::vector& visc_mult, std::vector& shear_mult) const; + + private: const PolymerProperties& polymer_props_; }; From 3bd59ab7e48923d81902c581528f97ca6c9cdc9a Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 10 Jun 2015 14:14:03 +0200 Subject: [PATCH 05/15] correcting the marco in the Point2D.hpp file. the typo is resulted from a stash apply CONFLICT. --- opm/polymer/Point2D.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/opm/polymer/Point2D.hpp b/opm/polymer/Point2D.hpp index 49549bef5..8d90b0990 100644 --- a/opm/polymer/Point2D.hpp +++ b/opm/polymer/Point2D.hpp @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef OPM_POINT2DINTERSECTION_INCLUDED -#define OPM_POINT2DINTERSECTION_INCLUDED +#ifndef OPM_POINT2D_HEADER_INCLUDED +#define OPM_POINT2D_HEADER_INCLUDED namespace Opm { @@ -98,5 +98,5 @@ namespace Opm { } // namespace Opm -#endif // OPM_POINT2DINTERSECTION_INCLUDED +#endif // OPM_POINT2D_HEADER_INCLUDED From df3c6cb8d1c53f78345f09355e644aaee8fc8e38 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 10 Jun 2015 14:40:46 +0200 Subject: [PATCH 06/15] using Water Position instead of assuming it is 0. --- opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 80ee16d67..b157a50cc 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -752,9 +752,10 @@ namespace Opm { // applying the shear-thinning to the water face if (has_plyshlog_) { + const int water_pos = pu.phase_pos[Water]; V shear_mult_wells_v = Eigen::Map(shear_mult_wells_.data(), shear_mult_wells_.size()); ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v); - mob_perfcells[0] = mob_perfcells[0] / shear_mult_wells_adb; + mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb; } @@ -901,7 +902,7 @@ namespace Opm { std::vector b_faces; - const int phase = fluid_.phaseUsage().phase_pos[BlackoilPhases::Aqua]; // water position + const int phase = fluid_.phaseUsage().phase_pos[Water]; // water position const int canonicalPhaseIdx = canph_[phase]; @@ -1100,8 +1101,9 @@ namespace Opm { } water_vel_wells.resize(cq_s[0].size()); - std::copy(&(cq_s[0].value()[0]), &(cq_s[0].value()[0]) + cq_s[0].size(), water_vel_wells.begin()); + const int water_pos = pu.phase_pos[Water]; + std::copy(&(cq_s[water_pos].value()[0]), &(cq_s[water_pos].value()[0]) + cq_s[water_pos].size(), water_vel_wells.begin()); const V& polymer_conc = state.concentration.value(); From 2b6163b5ff3642b777e47b995f9b52d67de98b30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 17 Jun 2015 12:51:39 +0200 Subject: [PATCH 07/15] Use the base class model_param_ and avoid recreating it. --- .../SimulatorFullyImplicitBlackoilPolymer_impl.hpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp index 4957ee6cf..94a75bc8e 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp @@ -64,12 +64,9 @@ namespace Opm -> std::unique_ptr { typedef typename Traits::Model Model; - typedef typename Model::ModelParameters ModelParams; - ModelParams modelParams( BaseType::param_ ); - typedef NewtonSolver Solver; - auto model = std::unique_ptr(new Model(modelParams, + auto model = std::unique_ptr(new Model(BaseType::model_param_, BaseType::grid_, BaseType::props_, BaseType::geo_, From 32590ec07294acc012481013f365c35a7018d46a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 18 Jun 2015 11:09:57 +0200 Subject: [PATCH 08/15] Use solver_param_ from base class. --- .../SimulatorFullyImplicitBlackoilPolymer_impl.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp index 94a75bc8e..8387819b6 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp @@ -86,9 +86,7 @@ namespace Opm model->setThresholdPressures(BaseType::threshold_pressures_by_face_); } - typedef typename Solver::SolverParameters SolverParams; - SolverParams solverParams( BaseType::param_ ); - return std::unique_ptr(new Solver(solverParams, std::move(model))); + return std::unique_ptr(new Solver(BaseType::solver_param_, std::move(model))); } From 2319420f06460b7cb4767ee3b61a026c6d910946 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 18 Jun 2015 14:34:52 +0200 Subject: [PATCH 09/15] Move polymer code extraAddWellEq -> addWellContributionToMassBalanceEq. --- .../fullyimplicit/BlackoilPolymerModel.hpp | 9 +++----- .../BlackoilPolymerModel_impl.hpp | 22 ++++++++++--------- 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index ccdc3a04f..09fbb214c 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -207,12 +207,9 @@ namespace Opm { assembleMassBalanceEq(const SolutionState& state); void - extraAddWellEq(const SolutionState& state, - const WellState& xw, - const std::vector& cq_ps, - const std::vector& cmix_s, - const ADB& cqt_is, - const std::vector& well_cells); + addWellContributionToMassBalanceEq(const SolutionState& state, + const WellState& xw, + const std::vector& cq_s); void computeMassFlux(const int actph , diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 71d97a68c..f3d8e2c0c 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -279,23 +279,25 @@ namespace Opm { template - void BlackoilPolymerModel::extraAddWellEq(const SolutionState& state, - const WellState& xw, - const std::vector& cq_ps, - const std::vector& cmix_s, - const ADB& cqt_is, - const std::vector& well_cells) + void BlackoilPolymerModel::addWellContributionToMassBalanceEq(const SolutionState& state, + const WellState& xw, + const std::vector& cq_s) { + Base::addWellContributionToMassBalanceEq(state, xw, cq_s); + // Add well contributions to polymer mass balance equation if (has_polymer_) { const ADB mc = computeMc(state); const int nc = xw.polymerInflow().size(); const V polyin = Eigen::Map(xw.polymerInflow().data(), nc); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); const V poly_in_perf = subset(polyin, well_cells); - const V poly_mc_perf = subset(mc, well_cells).value(); - const PhaseUsage& pu = fluid_.phaseUsage(); - const ADB cq_s_poly = cq_ps[pu.phase_pos[Water]] * poly_mc_perf - + cmix_s[pu.phase_pos[Water]] * cqt_is * poly_in_perf; + const V poly_mc_perf = subset(mc.value(), well_cells); + const ADB& cq_s_water = cq_s[fluid_.phaseUsage().phase_pos[Water]]; + Selector injector_selector(cq_s_water.value()); + const V poly_perf = injector_selector.select(poly_in_perf, poly_mc_perf); + const ADB cq_s_poly = cq_s_water * poly_perf; residual_.material_balance_eq[poly_pos_] -= superset(cq_s_poly, well_cells, nc); } } From 60494ac5310e4d181fa86b579faac4de142d072f Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 18 Jun 2015 16:56:39 +0200 Subject: [PATCH 10/15] changing to adapt to the change in the autodiff. --- .../fullyimplicit/BlackoilPolymerModel.hpp | 6 +- .../BlackoilPolymerModel_impl.hpp | 59 +++++++++++-------- 2 files changed, 40 insertions(+), 25 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 40971911b..ae16afcde 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -229,7 +229,11 @@ namespace Opm { void addWellEq(const SolutionState& state, WellState& xw, - V& aliveWells); + std::vector& mob_perfcells, + const std::vector& b_perfcells, + V& aliveWells, + std::vector& cq_s, + const bool welleq_initial); void addWellContributionToMassBalanceEq(const SolutionState& state, diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 2b59078ab..b8ab722ba 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -571,8 +571,30 @@ namespace Opm { assembleMassBalanceEq(state); // -------- Well equations ---------- + if ( ! wellsActive() ) { + return; + } + V aliveWells; + const int np = wells().number_of_phases; + std::vector cq_s(np, ADB::null()); + + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); + + std::vector mob_perfcells(np, ADB::null()); + std::vector b_perfcells(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + mob_perfcells[phase] = subset(rq_[phase].mob, well_cells); + b_perfcells[phase] = subset(rq_[phase].b, well_cells); + } + if (param_.solve_welleq_initially_ && initial_assembly) { + // solve the well equations as a pre-processing step + Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); + } + if (has_plyshlog_) { std::vector water_vel_wells; std::vector visc_mult_wells; @@ -594,7 +616,9 @@ namespace Opm { mob_perfcells = subset(rq_[0].mob,well_cells); */ } - addWellEq(state, well_state, aliveWells); + // addWellEq(state, well_state, aliveWells); + addWellEq(state, well_state, mob_perfcells, b_perfcells, aliveWells, cq_s, false); + addWellContributionToMassBalanceEq(state, well_state, cq_s); addWellControlEq(state, well_state, aliveWells); } @@ -725,12 +749,15 @@ namespace Opm { template void BlackoilPolymerModel::addWellEq(const SolutionState& state, - WellState& xw, - V& aliveWells) + WellState& xw, + std::vector& mob_perfcells, + const std::vector& b_perfcells, + V& aliveWells, + std::vector& cq_s, + const bool welleq_initial) { if( ! wellsActive() ) return ; - const int nc = Opm::AutoDiffGrid::numCells(grid_); const int np = wells().number_of_phases; const int nw = wells().number_of_wells; const int nperf = wells().well_connpos[nw]; @@ -740,27 +767,19 @@ namespace Opm { // pressure diffs computed already (once per step, not changing per iteration) const V& cdp = well_perforation_pressure_diffs_; - // Extract needed quantities for the perforation cells const ADB& p_perfcells = subset(state.pressure, well_cells); - const ADB& rv_perfcells = subset(state.rv,well_cells); - const ADB& rs_perfcells = subset(state.rs,well_cells); - std::vector mob_perfcells(np, ADB::null()); - std::vector b_perfcells(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - mob_perfcells[phase] = subset(rq_[phase].mob,well_cells); - b_perfcells[phase] = subset(rq_[phase].b,well_cells); - } + const ADB& rv_perfcells = subset(state.rv, well_cells); + const ADB& rs_perfcells = subset(state.rs, well_cells); // applying the shear-thinning to the water face - if (has_plyshlog_) { + if (has_plyshlog_ && (!welleq_initial)) { const int water_pos = pu.phase_pos[Water]; V shear_mult_wells_v = Eigen::Map(shear_mult_wells_.data(), shear_mult_wells_.size()); ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v); mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb; } - // Perforation pressure const ADB perfpressure = (wops_.w2p * state.bhp) + cdp; std::vector perfpressure_d(perfpressure.value().data(), perfpressure.value().data() + nperf); @@ -784,7 +803,6 @@ namespace Opm { } // HANDLE FLOW INTO WELLBORE - // compute phase volumetric rates at standard conditions std::vector cq_ps(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { @@ -801,7 +819,6 @@ namespace Opm { } // HANDLE FLOW OUT FROM WELLBORE - // Using total mobilities ADB total_mob = mob_perfcells[0]; for (int phase = 1; phase < np; ++phase) { @@ -826,7 +843,6 @@ namespace Opm { wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(V::Zero(nw)))) - q_ps; wbqt += wbq[phase]; } - // compute wellbore mixture at standard conditions. Selector notDeadWells_selector(wbqt.value(), Selector::Zero); std::vector cmix_s(np, ADB::null()); @@ -856,16 +872,10 @@ namespace Opm { ADB cqt_is = cqt_i/volumeRatio; // connection phase volumerates at standard conditions - std::vector cq_s(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; } - // Add well contributions to mass balance equations - for (int phase = 0; phase < np; ++phase) { - residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc); - } - // WELL EQUATIONS ADB qs = state.qs; for (int phase = 0; phase < np; ++phase) { @@ -893,6 +903,7 @@ namespace Opm { residual_.well_flux_eq = qs; } + template void BlackoilPolymerModel::computeWaterShearVelocityFaces(const V& transi, const std::vector& kr, From 4807a90c35ada08ded7747473206268433dc7cb0 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 19 Jun 2015 15:17:59 +0200 Subject: [PATCH 11/15] moving the shear-thinning application out of the addWellEq() So we still use the addWellEq() from the Base Class, without making a new addWellEq() function in the BlackoilPolymerModel class. --- .../fullyimplicit/BlackoilPolymerModel.hpp | 11 +- .../BlackoilPolymerModel_impl.hpp | 167 +----------------- 2 files changed, 8 insertions(+), 170 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index ae16afcde..a69df50bb 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -197,7 +197,7 @@ namespace Opm { using Base::drMaxRel; using Base::maxResidualAllowed; - // using Base::addWellEq; + using Base::addWellEq; using Base::updateWellControls; using Base::computeWellConnectionPressures; using Base::addWellControlEq; @@ -226,15 +226,6 @@ namespace Opm { void assembleMassBalanceEq(const SolutionState& state); - void - addWellEq(const SolutionState& state, - WellState& xw, - std::vector& mob_perfcells, - const std::vector& b_perfcells, - V& aliveWells, - std::vector& cq_s, - const bool welleq_initial); - void addWellContributionToMassBalanceEq(const SolutionState& state, const WellState& xw, diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index b8ab722ba..19bb9be38 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -606,6 +606,12 @@ namespace Opm { OPM_THROW(std::runtime_error, " failed in calculating the shear factors for wells "); } + // applying the shear-thinning to the water face + const int water_pos = fluid_.phaseUsage().phase_pos[Water]; + V shear_mult_wells_v = Eigen::Map(shear_mult_wells_.data(), shear_mult_wells_.size()); + ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v); + mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb; + /* const int nw = wells().number_of_wells; const int nperf = wells().well_connpos[nw]; @@ -617,7 +623,7 @@ namespace Opm { } // addWellEq(state, well_state, aliveWells); - addWellEq(state, well_state, mob_perfcells, b_perfcells, aliveWells, cq_s, false); + addWellEq(state, well_state, mob_perfcells, b_perfcells, aliveWells, cq_s); addWellContributionToMassBalanceEq(state, well_state, cq_s); addWellControlEq(state, well_state, aliveWells); } @@ -745,165 +751,6 @@ namespace Opm { return polymer_props_ad_.polymerWaterVelocityRatio(state.concentration); } - - template - void - BlackoilPolymerModel::addWellEq(const SolutionState& state, - WellState& xw, - std::vector& mob_perfcells, - const std::vector& b_perfcells, - V& aliveWells, - std::vector& cq_s, - const bool welleq_initial) - { - if( ! wellsActive() ) return ; - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - V Tw = Eigen::Map(wells().WI, nperf); - const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); - - // pressure diffs computed already (once per step, not changing per iteration) - const V& cdp = well_perforation_pressure_diffs_; - // Extract needed quantities for the perforation cells - const ADB& p_perfcells = subset(state.pressure, well_cells); - const ADB& rv_perfcells = subset(state.rv, well_cells); - const ADB& rs_perfcells = subset(state.rs, well_cells); - - // applying the shear-thinning to the water face - if (has_plyshlog_ && (!welleq_initial)) { - const int water_pos = pu.phase_pos[Water]; - V shear_mult_wells_v = Eigen::Map(shear_mult_wells_.data(), shear_mult_wells_.size()); - ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v); - mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb; - } - - // Perforation pressure - const ADB perfpressure = (wops_.w2p * state.bhp) + cdp; - std::vector perfpressure_d(perfpressure.value().data(), perfpressure.value().data() + nperf); - xw.perfPress() = perfpressure_d; - - // Pressure drawdown (also used to determine direction of flow) - const ADB drawdown = p_perfcells - perfpressure; - - // Compute vectors with zero and ones that - // selects the wanted quantities. - - // selects injection perforations - V selectInjectingPerforations = V::Zero(nperf); - // selects producing perforations - V selectProducingPerforations = V::Zero(nperf); - for (int c = 0; c < nperf; ++c){ - if (drawdown.value()[c] < 0) - selectInjectingPerforations[c] = 1; - else - selectProducingPerforations[c] = 1; - } - - // HANDLE FLOW INTO WELLBORE - // compute phase volumetric rates at standard conditions - std::vector cq_ps(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); - cq_ps[phase] = b_perfcells[phase] * cq_p; - } - if (active_[Oil] && active_[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const ADB cq_psOil = cq_ps[oilpos]; - const ADB cq_psGas = cq_ps[gaspos]; - cq_ps[gaspos] += rs_perfcells * cq_psOil; - cq_ps[oilpos] += rv_perfcells * cq_psGas; - } - - // HANDLE FLOW OUT FROM WELLBORE - // Using total mobilities - ADB total_mob = mob_perfcells[0]; - for (int phase = 1; phase < np; ++phase) { - total_mob += mob_perfcells[phase]; - } - // injection perforations total volume rates - const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); - - // compute wellbore mixture for injecting perforations - // The wellbore mixture depends on the inflow from the reservoar - // and the well injection rates. - - // compute avg. and total wellbore phase volumetric rates at standard conds - const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); - std::vector wbq(np, ADB::null()); - ADB wbqt = ADB::constant(V::Zero(nw)); - for (int phase = 0; phase < np; ++phase) { - const ADB& q_ps = wops_.p2w * cq_ps[phase]; - const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); - Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); - const int pos = pu.phase_pos[phase]; - wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(V::Zero(nw)))) - q_ps; - wbqt += wbq[phase]; - } - // compute wellbore mixture at standard conditions. - Selector notDeadWells_selector(wbqt.value(), Selector::Zero); - std::vector cmix_s(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const int pos = pu.phase_pos[phase]; - cmix_s[phase] = wops_.w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); - } - - // compute volume ratio between connection at standard conditions - ADB volumeRatio = ADB::constant(V::Zero(nperf)); - const ADB d = V::Constant(nperf,1.0) - rv_perfcells * rs_perfcells; - for (int phase = 0; phase < np; ++phase) { - ADB tmp = cmix_s[phase]; - - if (phase == Oil && active_[Gas]) { - const int gaspos = pu.phase_pos[Gas]; - tmp = tmp - rv_perfcells * cmix_s[gaspos] / d; - } - if (phase == Gas && active_[Oil]) { - const int oilpos = pu.phase_pos[Oil]; - tmp = tmp - rs_perfcells * cmix_s[oilpos] / d; - } - volumeRatio += tmp / b_perfcells[phase]; - } - - // injecting connections total volumerates at standard conditions - ADB cqt_is = cqt_i/volumeRatio; - - // connection phase volumerates at standard conditions - for (int phase = 0; phase < np; ++phase) { - cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; - } - - // WELL EQUATIONS - ADB qs = state.qs; - for (int phase = 0; phase < np; ++phase) { - qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); - - } - - // check for dead wells (used in the well controll equations) - aliveWells = V::Constant(nw, 1.0); - for (int w = 0; w < nw; ++w) { - if (wbqt.value()[w] == 0) { - aliveWells[w] = 0.0; - } - } - - // Update the perforation phase rates (used to calculate the pressure drop in the wellbore) - V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); - for (int phase = 1; phase < np; ++phase) { - cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); - } - - std::vector cq_d(cq.data(), cq.data() + nperf*np); - xw.perfPhaseRates() = cq_d; - - residual_.well_flux_eq = qs; - } - - template void BlackoilPolymerModel::computeWaterShearVelocityFaces(const V& transi, const std::vector& kr, From 3898ecde966b12434701c3bd58e7be69b6ce1db8 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Mon, 22 Jun 2015 15:43:01 +0800 Subject: [PATCH 12/15] Fix parameter errors. --- .../fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp index 293ed06f1..5823ff27a 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp @@ -29,6 +29,7 @@ #include #include #include +#include struct UnstructuredGrid; struct Wells; From 78775cd33e2abc8471cd813e96b26847b7ae781a Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 22 Jun 2015 12:48:22 +0200 Subject: [PATCH 13/15] parameter order change of addWellContributionToMassBalanceEq(). To follow the corresponding change in opm-autodiff. --- opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp | 6 +++--- opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp | 9 +++++---- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 09fbb214c..79507c843 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -207,9 +207,9 @@ namespace Opm { assembleMassBalanceEq(const SolutionState& state); void - addWellContributionToMassBalanceEq(const SolutionState& state, - const WellState& xw, - const std::vector& cq_s); + addWellContributionToMassBalanceEq(const std::vector& cq_s, + const SolutionState& state, + WellState& xw); void computeMassFlux(const int actph , diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index f3d8e2c0c..dbcaf149b 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -279,11 +279,12 @@ namespace Opm { template - void BlackoilPolymerModel::addWellContributionToMassBalanceEq(const SolutionState& state, - const WellState& xw, - const std::vector& cq_s) + void BlackoilPolymerModel::addWellContributionToMassBalanceEq(const std::vector& cq_s, + const SolutionState& state, + WellState& xw) + { - Base::addWellContributionToMassBalanceEq(state, xw, cq_s); + Base::addWellContributionToMassBalanceEq(cq_s, state, xw); // Add well contributions to polymer mass balance equation if (has_polymer_) { From 39986c3363ab3d1397e3917cc874883bff58d4b8 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 22 Jun 2015 15:38:31 +0200 Subject: [PATCH 14/15] adapting to the new refactoring in opm-autodiff#406. replacing addWellEq() with seperated functions, which make the incoropration of the shear-thinning much natural and easier. --- .../fullyimplicit/BlackoilPolymerModel.hpp | 6 +- .../BlackoilPolymerModel_impl.hpp | 194 ++++-------------- 2 files changed, 47 insertions(+), 153 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 1ce3f963c..c45caaa84 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -197,7 +197,6 @@ namespace Opm { using Base::drMaxRel; using Base::maxResidualAllowed; - using Base::addWellEq; using Base::updateWellControls; using Base::computeWellConnectionPressures; using Base::addWellControlEq; @@ -283,9 +282,8 @@ namespace Opm { const std::vector& phasePressure, const SolutionState& state, std::vector& water_vel, std::vector& visc_mult); - void computeWaterShearVelocityWells(const SolutionState& state, WellState& well_state, - V& aliveWells, std::vector& water_vel_wells, std::vector& visc_mult_wells); - + void computeWaterShearVelocityWells(const SolutionState& state, WellState& xw, const ADB& cq_sw, + std::vector& water_vel_wells, std::vector& visc_mult_wells); }; diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index e4e2480b1..5d22916d7 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -287,7 +287,7 @@ namespace Opm { std::vector visc_mult; computeWaterShearVelocityFaces(transi, kr, state.canonical_phase_pressures, state, water_vel, visc_mult); - if(!polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_)) { + if ( !polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_) ) { // std::cerr << " failed in calculating the shear-multiplier " << std::endl; OPM_THROW(std::runtime_error, " failed in calculating the shear-multiplier. "); } @@ -596,19 +596,21 @@ namespace Opm { Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); } + Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); + if (has_plyshlog_) { std::vector water_vel_wells; std::vector visc_mult_wells; - computeWaterShearVelocityWells(state, well_state, aliveWells, water_vel_wells, visc_mult_wells); + const int water_pos = fluid_.phaseUsage().phase_pos[Water]; + computeWaterShearVelocityWells(state, well_state, cq_s[water_pos], water_vel_wells, visc_mult_wells); - if (!polymer_props_ad_.computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_)) { + if ( !polymer_props_ad_.computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_) ) { // std::cout << " failed in calculating the shear factors for wells " << std::endl; OPM_THROW(std::runtime_error, " failed in calculating the shear factors for wells "); } - // applying the shear-thinning to the water face - const int water_pos = fluid_.phaseUsage().phase_pos[Water]; + // applying the shear-thinning to the water phase V shear_mult_wells_v = Eigen::Map(shear_mult_wells_.data(), shear_mult_wells_.size()); ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v); mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb; @@ -623,9 +625,10 @@ namespace Opm { mob_perfcells = subset(rq_[0].mob,well_cells); */ } - // addWellEq(state, well_state, aliveWells); - addWellEq(state, well_state, mob_perfcells, b_perfcells, aliveWells, cq_s); - addWellContributionToMassBalanceEq(state, well_state, cq_s); + Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); + Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state); + Base::addWellFluxEq(cq_s, state); + addWellContributionToMassBalanceEq(cq_s, state, well_state); addWellControlEq(state, well_state, aliveWells); } @@ -802,9 +805,9 @@ namespace Opm { rq_[ phase ].mflux = upwind.select(b * mob) * (transi * dh); - const auto& tempb_faces = upwind.select(b); - b_faces.resize(tempb_faces.size()); - std::copy(&(tempb_faces.value()[0]), &(tempb_faces.value()[0]) + tempb_faces.size(), b_faces.begin()); + const auto& b_faces_adb = upwind.select(b); + b_faces.resize(b_faces_adb.size()); + std::copy(&(b_faces_adb.value()[0]), &(b_faces_adb.value()[0]) + b_faces_adb.size(), b_faces.begin()); const auto& internal_faces = ops_.internal_faces; @@ -816,11 +819,11 @@ namespace Opm { } const ADB phi = Opm::AutoDiffBlock::constant(Eigen::Map(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1)); - const ADB temp_phiavg = ops_.caver * phi; + const ADB phiavg_adb = ops_.caver * phi; std::vector phiavg; - phiavg.resize(temp_phiavg.size()); - std::copy(&(temp_phiavg.value()[0]), &(temp_phiavg.value()[0]) + temp_phiavg.size(), phiavg.begin()); + phiavg.resize(phiavg_adb.size()); + std::copy(&(phiavg_adb.value()[0]), &(phiavg_adb.value()[0]) + phiavg_adb.size(), phiavg.begin()); water_vel.resize(nface); std::copy(&(rq_[0].mflux.value()[0]), &(rq_[0].mflux.value()[0]) + nface, water_vel.begin()); @@ -832,146 +835,42 @@ namespace Opm { template void - BlackoilPolymerModel::computeWaterShearVelocityWells(const SolutionState& state, WellState& xw, - V& aliveWells, std::vector& water_vel_wells, std::vector& visc_mult_wells) + BlackoilPolymerModel::computeWaterShearVelocityWells(const SolutionState& state, WellState& xw, const ADB& cq_sw, + std::vector& water_vel_wells, std::vector& visc_mult_wells) { if( ! wellsActive() ) return ; - // const int nc = Opm::AutoDiffGrid::numCells(grid_); - const int np = wells().number_of_phases; const int nw = wells().number_of_wells; const int nperf = wells().well_connpos[nw]; - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - V Tw = Eigen::Map(wells().WI, nperf); const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); - // pressure diffs computed already (once per step, not changing per iteration) - const V& cdp = well_perforation_pressure_diffs_; - - // Extract needed quantities for the perforation cells - const ADB& p_perfcells = subset(state.pressure, well_cells); - const ADB& rv_perfcells = subset(state.rv,well_cells); - const ADB& rs_perfcells = subset(state.rs,well_cells); - std::vector mob_perfcells(np, ADB::null()); - std::vector b_perfcells(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - mob_perfcells[phase] = subset(rq_[phase].mob,well_cells); - b_perfcells[phase] = subset(rq_[phase].b,well_cells); - } - - // Perforation pressure - const ADB perfpressure = (wops_.w2p * state.bhp) + cdp; - std::vector perfpressure_d(perfpressure.value().data(), perfpressure.value().data() + nperf); - xw.perfPress() = perfpressure_d; - - // Pressure drawdown (also used to determine direction of flow) - const ADB drawdown = p_perfcells - perfpressure; - - // Compute vectors with zero and ones that - // selects the wanted quantities. - - // selects injection perforations - V selectInjectingPerforations = V::Zero(nperf); - // selects producing perforations - V selectProducingPerforations = V::Zero(nperf); - for (int c = 0; c < nperf; ++c) { - if (drawdown.value()[c] < 0) - selectInjectingPerforations[c] = 1; - else - selectProducingPerforations[c] = 1; - } - - // HANDLE FLOW INTO WELLBORE - - // compute phase volumetric rates at standard conditions - std::vector cq_ps(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); - cq_ps[phase] = b_perfcells[phase] * cq_p; - } - if (active_[Oil] && active_[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const ADB cq_psOil = cq_ps[oilpos]; - const ADB cq_psGas = cq_ps[gaspos]; - cq_ps[gaspos] += rs_perfcells * cq_psOil; - cq_ps[oilpos] += rv_perfcells * cq_psGas; - } - - // HANDLE FLOW OUT FROM WELLBORE - - // Using total mobilities - ADB total_mob = mob_perfcells[0]; - for (int phase = 1; phase < np; ++phase) { - total_mob += mob_perfcells[phase]; - } - // injection perforations total volume rates - const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); - - // compute wellbore mixture for injecting perforations - // The wellbore mixture depends on the inflow from the reservoar - // and the well injection rates. - - // compute avg. and total wellbore phase volumetric rates at standard conds - const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); - std::vector wbq(np, ADB::null()); - ADB wbqt = ADB::constant(V::Zero(nw)); - for (int phase = 0; phase < np; ++phase) { - const ADB& q_ps = wops_.p2w * cq_ps[phase]; - const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); - Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); - const int pos = pu.phase_pos[phase]; - wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(V::Zero(nw)))) - q_ps; - wbqt += wbq[phase]; - } - - // compute wellbore mixture at standard conditions. - Selector notDeadWells_selector(wbqt.value(), Selector::Zero); - std::vector cmix_s(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const int pos = pu.phase_pos[phase]; - cmix_s[phase] = wops_.w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); - } - - // compute volume ratio between connection at standard conditions - ADB volumeRatio = ADB::constant(V::Zero(nperf)); - const ADB d = V::Constant(nperf,1.0) - rv_perfcells * rs_perfcells; - for (int phase = 0; phase < np; ++phase) { - ADB tmp = cmix_s[phase]; - - if (phase == Oil && active_[Gas]) { - const int gaspos = pu.phase_pos[Gas]; - tmp = tmp - rv_perfcells * cmix_s[gaspos] / d; - } - if (phase == Gas && active_[Oil]) { - const int oilpos = pu.phase_pos[Oil]; - tmp = tmp - rs_perfcells * cmix_s[oilpos] / d; - } - volumeRatio += tmp / b_perfcells[phase]; - } - - // injecting connections total volumerates at standard conditions - ADB cqt_is = cqt_i/volumeRatio; - - // connection phase volumerates at standard conditions - std::vector cq_s(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; - } - - water_vel_wells.resize(cq_s[0].size()); - - const int water_pos = pu.phase_pos[Water]; - std::copy(&(cq_s[water_pos].value()[0]), &(cq_s[water_pos].value()[0]) + cq_s[water_pos].size(), water_vel_wells.begin()); + water_vel_wells.resize(cq_sw.size()); + std::copy(&(cq_sw.value()[0]), &(cq_sw.value()[0]) + cq_sw.size(), water_vel_wells.begin()); const V& polymer_conc = state.concentration.value(); V visc_mult_cells = polymer_props_ad_.viscMult(polymer_conc); + V visc_mult_wells_v = subset(visc_mult_cells, well_cells); - V temp_visc_mult_wells = subset(visc_mult_cells, well_cells); + visc_mult_wells.resize(visc_mult_wells_v.size()); + std::copy(&(visc_mult_wells_v[0]), &(visc_mult_wells_v[0]) + visc_mult_wells_v.size(), visc_mult_wells.begin()); - visc_mult_wells.resize(temp_visc_mult_wells.size()); - std::copy(&(temp_visc_mult_wells[0]), &(temp_visc_mult_wells[0]) + temp_visc_mult_wells.size(), visc_mult_wells.begin()); + const int water_pos = fluid_.phaseUsage().phase_pos[Water]; + ADB b_perfcells = subset(rq_[water_pos].b, well_cells); + + const ADB& p_perfcells = subset(state.pressure, well_cells); + const V& cdp = well_perforation_pressure_diffs_; + const ADB perfpressure = (wops_.w2p * state.bhp) + cdp; + // Pressure drawdown (also used to determine direction of flow) + const ADB drawdown = p_perfcells - perfpressure; + + // selects injection perforations + V selectInjectingPerforations = V::Zero(nperf); + for (int c = 0; c < nperf; ++c) { + if (drawdown.value()[c] < 0) { + selectInjectingPerforations[c] = 1; + } + } // for the injection wells for (int i = 0; i < well_cells.size(); ++i) { @@ -981,18 +880,15 @@ namespace Opm { } const ADB phi = Opm::AutoDiffBlock::constant(Eigen::Map(& fluid_.porosity()[0], AutoDiffGrid::numCells(grid_), 1)); - - const ADB temp_phi_wells = subset(phi, well_cells); + const ADB phi_wells_adb = subset(phi, well_cells); std::vector phi_wells; - phi_wells.resize(temp_phi_wells.size()); - - std::copy(&(temp_phi_wells.value()[0]), &(temp_phi_wells.value()[0]) + temp_phi_wells.size(), phi_wells.begin()); + phi_wells.resize(phi_wells_adb.size()); + std::copy(&(phi_wells_adb.value()[0]), &(phi_wells_adb.value()[0]) + phi_wells_adb.size(), phi_wells.begin()); std::vector b_wells; - - b_wells.resize(b_perfcells[0].size()); - std::copy(&(b_perfcells[0].value()[0]), &(b_perfcells[0].value()[0]) + b_perfcells[0].size(), b_wells.begin()); + b_wells.resize(b_perfcells.size()); + std::copy(&(b_perfcells.value()[0]), &(b_perfcells.value()[0]) + b_perfcells.size(), b_wells.begin()); for (int i = 0; i < water_vel_wells.size(); ++i) { water_vel_wells[i] = b_wells[i] * water_vel_wells[i] / (phi_wells[i] * 2. * M_PI * wells_rep_radius_[i] * wells_perf_length_[i]); From 5cdc6776727b63aaa7741dc06b33f47954ddc98d Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 22 Jun 2015 15:41:56 +0200 Subject: [PATCH 15/15] cleaning up commented codes and adding a few comments. --- opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp | 2 ++ .../fullyimplicit/BlackoilPolymerModel_impl.hpp | 10 ---------- 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index c45caaa84..5521800d6 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -282,6 +282,8 @@ namespace Opm { const std::vector& phasePressure, const SolutionState& state, std::vector& water_vel, std::vector& visc_mult); + /// Computing the water velocity without shear-thinning for the well perforations based on the water flux rate. + /// The water velocity will be used for shear-thinning calculation. void computeWaterShearVelocityWells(const SolutionState& state, WellState& xw, const ADB& cq_sw, std::vector& water_vel_wells, std::vector& visc_mult_wells); diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 5d22916d7..9614e61d2 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -606,7 +606,6 @@ namespace Opm { computeWaterShearVelocityWells(state, well_state, cq_s[water_pos], water_vel_wells, visc_mult_wells); if ( !polymer_props_ad_.computeShearMultLog(water_vel_wells, visc_mult_wells, shear_mult_wells_) ) { - // std::cout << " failed in calculating the shear factors for wells " << std::endl; OPM_THROW(std::runtime_error, " failed in calculating the shear factors for wells "); } @@ -614,15 +613,6 @@ namespace Opm { V shear_mult_wells_v = Eigen::Map(shear_mult_wells_.data(), shear_mult_wells_.size()); ADB shear_mult_wells_adb = ADB::constant(shear_mult_wells_v); mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb; - - /* const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - - const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); - - // assuming the water phase is the first phase - const int nw = wells().number_of_wells; - mob_perfcells = subset(rq_[0].mob,well_cells); */ } Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s);