From 97f5c5ace5585780873552122763bf6fd86a72f0 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Fri, 10 Jan 2014 17:20:08 +0800 Subject: [PATCH] fix bug when computing the maxconcentration that the cell experienced. --- ...FullyImplicitCompressiblePolymerSolver.cpp | 27 ++++++++++--------- ...FullyImplicitCompressiblePolymerSolver.hpp | 4 +-- .../FullyImplicitTwophasePolymerSolver.cpp | 25 +++++++++-------- .../FullyImplicitTwophasePolymerSolver.hpp | 8 +++--- 4 files changed, 31 insertions(+), 33 deletions(-) diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 96ce6bcdc..e7f71dbb6 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -297,7 +297,7 @@ namespace { FullyImplicitCompressiblePolymerSolver::SolutionState FullyImplicitCompressiblePolymerSolver::constantState(const PolymerBlackoilState& x, - const WellState& xw) + const WellState& xw) { const int nc = grid_.number_of_cells; const int np = x.numPhases(); @@ -342,7 +342,7 @@ namespace { const int nw = wells_.number_of_wells; // The transpose() below switches the ordering. const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); - const V qs = Eigen::Map(wrates.data(), nw*np); + const V qs = Eigen::Map(wrates.data(), nw * np); state.qs = ADB::constant(qs, bpat); // Well bottom-hole pressure. @@ -359,7 +359,7 @@ namespace { FullyImplicitCompressiblePolymerSolver::SolutionState FullyImplicitCompressiblePolymerSolver::variableState(const PolymerBlackoilState& x, - const WellState& xw) + const WellState& xw) { const int nc = grid_.number_of_cells; const int np = x.numPhases(); @@ -435,7 +435,7 @@ namespace { void FullyImplicitCompressiblePolymerSolver::computeAccum(const SolutionState& state, - const int aix ) + const int aix ) { const ADB& press = state.pressure; @@ -454,18 +454,18 @@ namespace { - ADB + V FullyImplicitCompressiblePolymerSolver:: - computeCmax(const ADB& c) const + computeCmax(const PolymerBlackoilState& x) const { const int nc = c.value().size(); - V cmax(nc); - - for (int i = 0; i < nc; ++i) { - cmax(i) = (cmax(i) > c.value()(i)) ? cmax(i) : c.value()(i); + const V cmax = Eigen::Map(& x.maxconcentration()[0], nc); + const V c = Eigen::Map(& x.concentration()[0], nc); + for (int i = 0; i < nc; ++i) { + cmax(i) = std::max(cmax(i), c(i)); } - - return ADB::constant(cmax, c.blockPattern()); + + return cmax; } void @@ -493,7 +493,8 @@ namespace { const V trans = subset(geo_.transmissibility(), ops_.internal_faces); const std::vector kr = computeRelPerm(state); const double dead_pore_vol = polymer_props_ad_.deadPoreVol(); - const ADB cmax = computeCmax(state.concentration); + const V cmax_v = computeCmax(x); + const ADB cmax = ADB::constant(cmax_v, state.concentration.blockPattern()); const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax); const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]); const ADB mc = computeMc(state); diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp index 27f9596c6..77d1348ec 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp @@ -191,8 +191,8 @@ namespace Opm { computeFracFlow(const ADB& kro, const ADB& krw_eff, const ADB& c) const; - ADB - computeCmax(const ADB& c) const; + V + computeCmax(const PolymerBlackoilState& x) const; ADB computeMc(const SolutionState& state) const; ADB diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp index 98a78deee..a1e4fa730 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp @@ -1,4 +1,3 @@ - #include #include @@ -304,7 +303,7 @@ namespace { const int nw = wells_.number_of_wells; // The transpose() below switches the ordering. const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); - const V qs = Eigen::Map(wrates.data(), nw*np); + const V qs = Eigen::Map(wrates.data(), nw * np); vars0.push_back(qs); // Initial well bottom hole pressure. @@ -335,6 +334,7 @@ namespace { // Qs. state.qs = vars[ nextvar++ ]; + // BHP. state.bhp = vars[ nextvar++ ]; assert(nextvar == int(vars.size())); @@ -343,18 +343,18 @@ namespace { } - ADB + V FullyImplicitTwophasePolymerSolver:: - computeCmax(const ADB& c) const + computeCmax(const PolymerState& x) const { const int nc = c.value().size(); - V cmax(nc); - + const V cmax = Eigen::Map(& x.maxconcentration()[0], nc); + const V c = Eigen::Map(& x.concentration()[0], nc); for (int i = 0; i < nc; ++i) { - cmax(i) = (cmax(i) > c.value()(i)) ? cmax(i) : c.value()(i); + cmax(i) = std::max(cmax(i), c(i)); } - return ADB::constant(cmax, c.blockPattern()); + return cmax; } @@ -374,7 +374,8 @@ namespace { const V trans = subset(transmissibility(), ops_.internal_faces); const std::vector kr = computeRelPerm(state); - const ADB cmax = computeCmax(state.concentration); + const V cmax_v = computeCmax(x); + const ADB cmax = ADB::constant(cmax, state.concentration.blockPattern()); const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax); const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]); const ADB mc = computeMc(state); @@ -594,7 +595,7 @@ namespace { const V inSrc = Eigen::Map(& insrc[0], grid_.number_of_cells); const V polyin = Eigen::Map(& polymer_inflow_c[0], grid_.number_of_cells); // compute the out-fracflow. - const std::vector f = computeFracFlow(kro, krw_eff, c); + const std::vector f = computeFracFlow(); // compute the in-fracflow. V zero = V::Zero(grid_.number_of_cells); V one = V::Ones(grid_.number_of_cells); @@ -614,9 +615,7 @@ namespace { std::vector - FullyImplicitTwophasePolymerSolver::computeFracFlow(const ADB& kro, - const ADB& krw_eff, - const ADB& c) const + FullyImplicitTwophasePolymerSolver::computeFracFlow() const { ADB total_mob = mob_[0] + mob_[1]; diff --git a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp index 865332ec1..1a6c5e438 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp @@ -105,9 +105,7 @@ namespace Opm { std::vector - computeFracFlow(const ADB& kro, - const ADB& krw_eff, - const ADB& c) const; + computeFracFlow() const; double residualNorm() const; ADB @@ -116,8 +114,8 @@ namespace Opm { const std::vector& polymer_inflow_c, const SolutionState& state) const; - ADB - computeCmax(const ADB& c) const; + V + computeCmax(const PolymerState& x) const; ADB computeMc(const SolutionState& state) const; ADB